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 561 by mmeineke, Fri Jun 20 20:29:36 2003 UTC vs.
Revision 563 by mmeineke, Mon Jun 23 21:24:03 2003 UTC

# Line 137 | Line 137 | void Integrator::checkConstraints( void ){
137        constrainedA[i] = temp_con[i].get_a();
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";
143      }
144  
145      
# Line 402 | Line 405 | void Integrator::preMove( void ){
405  
406    if( nConstrained ){
407  
405 – //    if( oldAtoms != nAtoms ){
406 –      
407 – //       // save oldAtoms to check for lode balanceing later on.
408 –      
409 – //       oldAtoms = nAtoms;
410 –      
411 – //       delete[] moving;
412 – //       delete[] moved;
413 – //       delete[] oldPos;
414 –      
415 – //       moving = new int[nAtoms];
416 – //       moved  = new int[nAtoms];
417 –      
418 – //       oldPos = new double[nAtoms*3];
419 – //     }
420 –  
408      for(i=0; i<(nAtoms*3); i++) oldPos[i] = pos[i];
409    }
410   }  
# Line 428 | Line 415 | void Integrator::constrainA(){
415    int done;
416    double pxab, pyab, pzab;
417    double rxab, ryab, rzab;
418 <  int a, b;
418 >  int a, b, ax, ay, az, bx, by, bz;
419    double rma, rmb;
420    double dx, dy, dz;
421    double rpab;
# Line 455 | Line 442 | void Integrator::constrainA(){
442  
443        a = constrainedA[i];
444        b = constrainedB[i];
445 <    
445 >      
446 >      ax = (a*3) + 0;
447 >      ay = (a*3) + 1;
448 >      az = (a*3) + 2;
449 >
450 >      bx = (b*3) + 0;
451 >      by = (b*3) + 1;
452 >      bz = (b*3) + 2;
453 >
454 >
455        if( moved[a] || moved[b] ){
456          
457 <        pxab = pos[3*a+0] - pos[3*b+0];
458 <        pyab = pos[3*a+1] - pos[3*b+1];
459 <        pzab = pos[3*a+2] - pos[3*b+2];
457 >        pxab = pos[ax] - pos[bx];
458 >        pyab = pos[ay] - pos[by];
459 >        pzab = pos[az] - pos[bz];
460  
461          //periodic boundary condition
462          pxab = pxab - info->box_x * copysign(1, pxab)
463 <          * int( fabs(pxab) / info->box_x + 0.5);
463 >          * (int)( fabs(pxab / info->box_x) + 0.5);
464          pyab = pyab - info->box_y * copysign(1, pyab)
465 <          * int( fabs(pyab) / info->box_y + 0.5);
465 >          * (int)( fabs(pyab / info->box_y) + 0.5);
466          pzab = pzab - info->box_z * copysign(1, pzab)
467 <          * int( fabs(pzab) / info->box_z + 0.5);
467 >          * (int)( fabs(pzab / info->box_z) + 0.5);
468        
469          pabsq = pxab * pxab + pyab * pyab + pzab * pzab;
470          rabsq = constrainedDsqr[i];
471          diffsq = pabsq - rabsq;
472  
473          // the original rattle code from alan tidesley
474 <        if (fabs(diffsq) > tol*rabsq*2) {
475 <          rxab = oldPos[3*a+0] - oldPos[3*b+0];
476 <          ryab = oldPos[3*a+1] - oldPos[3*b+1];
477 <          rzab = oldPos[3*a+2] - oldPos[3*b+2];
474 >        if (fabs(diffsq) > (tol*rabsq*2)) {
475 >          rxab = oldPos[ax] - oldPos[bx];
476 >          ryab = oldPos[ay] - oldPos[by];
477 >          rzab = oldPos[az] - oldPos[bz];
478  
479            rxab = rxab - info->box_x * copysign(1, rxab)
480 <            * int( fabs(rxab) / info->box_x + 0.5);
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);
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);
484 >            * (int)( fabs(rzab / info->box_z) + 0.5);
485  
486            rpab = rxab * pxab + ryab * pyab + rzab * pzab;
487            rpabsq = rpab * rpab;
488  
489  
490            if (rpabsq < (rabsq * -diffsq)){
491 +
492 +            cerr << "rpabsq = " << rpabsq << ", rabsq = " << rabsq
493 +                 << ", -diffsq = " << -diffsq << "\n";
494 +
495   #ifdef IS_MPI
496              a = atoms[a]->getGlobalIndex();
497              b = atoms[b]->getGlobalIndex();
498   #endif //is_mpi
499              sprintf( painCave.errMsg,
500 <                     "Constraint failure in constrainA at atom %d and %d\n.",
500 >                     "Constraint failure in constrainA at atom %d and %d.\n",
501                       a, b );
502              painCave.isFatal = 1;
503              simError();
# Line 511 | Line 511 | void Integrator::constrainA(){
511            dy = ryab * gab;
512            dz = rzab * gab;
513  
514 <          pos[3*a+0] += rma * dx;
515 <          pos[3*a+1] += rma * dy;
516 <          pos[3*a+2] += rma * dz;
514 >          pos[ax] += rma * dx;
515 >          pos[ay] += rma * dy;
516 >          pos[az] += rma * dz;
517  
518 <          pos[3*b+0] -= rmb * dx;
519 <          pos[3*b+1] -= rmb * dy;
520 <          pos[3*b+2] -= rmb * dz;
518 >          pos[bx] -= rmb * dx;
519 >          pos[by] -= rmb * dy;
520 >          pos[bz] -= rmb * dz;
521  
522            dx = dx / dt;
523            dy = dy / dt;
524            dz = dz / dt;
525  
526 <          vel[3*a+0] += rma * dx;
527 <          vel[3*a+1] += rma * dy;
528 <          vel[3*a+2] += rma * dz;
526 >          vel[ax] += rma * dx;
527 >          vel[ay] += rma * dy;
528 >          vel[az] += rma * dz;
529  
530 <          vel[3*b+0] -= rmb * dx;
531 <          vel[3*b+1] -= rmb * dy;
532 <          vel[3*b+2] -= rmb * dz;
530 >          vel[bx] -= rmb * dx;
531 >          vel[by] -= rmb * dy;
532 >          vel[bz] -= rmb * dz;
533  
534            moving[a] = 1;
535            moving[b] = 1;
# Line 545 | Line 545 | void Integrator::constrainA(){
545      }
546  
547      iteration++;
548 +    cerr << "iterainA = " << iteration << "\n";
549    }
550  
551    if( !done ){
# Line 564 | Line 565 | void Integrator::constrainB( void ){
565    int done;
566    double vxab, vyab, vzab;
567    double rxab, ryab, rzab;
568 <  int a, b;
568 >  int a, b, ax, ay, az, bx, by, bz;
569    double rma, rmb;
570    double dx, dy, dz;
571    double rabsq, pabsq, rvab;
# Line 586 | Line 587 | void Integrator::constrainB( void ){
587        a = constrainedA[i];
588        b = constrainedB[i];
589  
590 +      ax = 3*a +0;
591 +      ay = 3*a +1;
592 +      az = 3*a +2;
593 +
594 +      bx = 3*b +0;
595 +      by = 3*b +1;
596 +      bz = 3*b +2;
597 +
598        if( moved[a] || moved[b] ){
599          
600 <        vxab = vel[3*a+0] - vel[3*b+0];
601 <        vyab = vel[3*a+1] - vel[3*b+1];
602 <        vzab = vel[3*a+2] - vel[3*b+2];
600 >        vxab = vel[ax] - vel[bx];
601 >        vyab = vel[ay] - vel[by];
602 >        vzab = vel[az] - vel[bz];
603  
604 <        rxab = pos[3*a+0] - pos[3*b+0];
605 <        ryab = pos[3*a+1] - pos[3*b+1];
606 <        rzab = pos[3*a+2] - pos[3*b+2];
604 >        rxab = pos[ax] - pos[bx];
605 >        ryab = pos[ay] - pos[by];
606 >        rzab = pos[az] - pos[bz];
607          
608          rxab = rxab - info->box_x * copysign(1, rxab)
609 <          * int( fabs(rxab) / info->box_x + 0.5);
609 >          * (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);
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);
613 >          * (int)( fabs(rzab / info->box_z) + 0.5);
614  
615          rma = 1.0 / atoms[a]->getMass();
616          rmb = 1.0 / atoms[b]->getMass();
# Line 616 | Line 625 | void Integrator::constrainB( void ){
625            dy = ryab * gab;
626            dz = rzab * gab;
627            
628 <          vel[3*a+0] += rma * dx;
629 <          vel[3*a+1] += rma * dy;
630 <          vel[3*a+2] += rma * dz;
628 >          vel[ax] += rma * dx;
629 >          vel[ay] += rma * dy;
630 >          vel[az] += rma * dz;
631  
632 <          vel[3*b+0] -= rmb * dx;
633 <          vel[3*b+1] -= rmb * dy;
634 <          vel[3*b+2] -= rmb * dz;
632 >          vel[bx] -= rmb * dx;
633 >          vel[by] -= rmb * dy;
634 >          vel[bz] -= rmb * dz;
635            
636            moving[a] = 1;
637            moving[b] = 1;

Diff Legend

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