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 567 by mmeineke, Wed Jun 25 21:12:14 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      }
142  
143      
# Line 292 | Line 293 | void Integrator::moveA( void ){
293    double Tb[3];
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];
307 <
308 <  
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 402 | Line 403 | void Integrator::preMove( void ){
403  
404    if( nConstrained ){
405  
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 –  
406      for(i=0; i<(nAtoms*3); i++) oldPos[i] = pos[i];
407    }
408   }  
# Line 428 | Line 413 | void Integrator::constrainA(){
413    int done;
414    double pxab, pyab, pzab;
415    double rxab, ryab, rzab;
416 <  int a, b;
416 >  int a, b, ax, ay, az, bx, by, bz;
417    double rma, rmb;
418    double dx, dy, dz;
419    double rpab;
# Line 444 | Line 429 | void Integrator::constrainA(){
429      moving[i] = 0;
430      moved[i]  = 1;
431    }
432 <  
448 <  
432 >
433    iteration = 0;
434    done = 0;
435    while( !done && (iteration < maxIteration )){
# Line 455 | Line 439 | void Integrator::constrainA(){
439  
440        a = constrainedA[i];
441        b = constrainedB[i];
442 <    
442 >      
443 >      ax = (a*3) + 0;
444 >      ay = (a*3) + 1;
445 >      az = (a*3) + 2;
446 >
447 >      bx = (b*3) + 0;
448 >      by = (b*3) + 1;
449 >      bz = (b*3) + 2;
450 >
451        if( moved[a] || moved[b] ){
452          
453 <        pxab = pos[3*a+0] - pos[3*b+0];
454 <        pyab = pos[3*a+1] - pos[3*b+1];
455 <        pzab = pos[3*a+2] - pos[3*b+2];
453 >        pxab = pos[ax] - pos[bx];
454 >        pyab = pos[ay] - pos[by];
455 >        pzab = pos[az] - pos[bz];
456  
457 <        //periodic boundary condition
457 >        //periodic boundary condition
458          pxab = pxab - info->box_x * copysign(1, pxab)
459 <          * int( fabs(pxab) / info->box_x + 0.5);
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);
461 >          * (int)( fabs(pyab / info->box_y) + 0.5);
462          pzab = pzab - info->box_z * copysign(1, pzab)
463 <          * int( fabs(pzab) / info->box_z + 0.5);
464 <      
465 <        pabsq = pxab * pxab + pyab * pyab + pzab * pzab;
463 >          * (int)( fabs(pzab / info->box_z) + 0.5);
464 >
465 >        pabsq = pxab * pxab + pyab * pyab + pzab * pzab;
466 >
467          rabsq = constrainedDsqr[i];
468 <        diffsq = pabsq - rabsq;
468 >        diffsq = rabsq - pabsq;
469  
470          // the original rattle code from alan tidesley
471 <        if (fabs(diffsq) > tol*rabsq*2) {
472 <          rxab = oldPos[3*a+0] - oldPos[3*b+0];
473 <          ryab = oldPos[3*a+1] - oldPos[3*b+1];
474 <          rzab = oldPos[3*a+2] - oldPos[3*b+2];
475 <
471 >        if (fabs(diffsq) > (tol*rabsq*2)) {
472 >          rxab = oldPos[ax] - oldPos[bx];
473 >          ryab = oldPos[ay] - oldPos[by];
474 >          rzab = oldPos[az] - oldPos[bz];
475 >
476            rxab = rxab - info->box_x * copysign(1, rxab)
477 <            * int( fabs(rxab) / info->box_x + 0.5);
477 >            * (int)( fabs(rxab / info->box_x) + 0.5);
478            ryab = ryab - info->box_y * copysign(1, ryab)
479 <            * int( fabs(ryab) / info->box_y + 0.5);
479 >            * (int)( fabs(ryab / info->box_y) + 0.5);
480            rzab = rzab - info->box_z * copysign(1, rzab)
481 <            * int( fabs(rzab) / info->box_z + 0.5);
481 >            * (int)( fabs(rzab / info->box_z) + 0.5);
482  
483            rpab = rxab * pxab + ryab * pyab + rzab * pzab;
484 +
485            rpabsq = rpab * rpab;
486  
487  
488            if (rpabsq < (rabsq * -diffsq)){
489 +
490   #ifdef IS_MPI
491              a = atoms[a]->getGlobalIndex();
492              b = atoms[b]->getGlobalIndex();
493   #endif //is_mpi
494              sprintf( painCave.errMsg,
495 <                     "Constraint failure in constrainA at atom %d and %d\n.",
495 >                     "Constraint failure in constrainA at atom %d and %d.\n",
496                       a, b );
497              painCave.isFatal = 1;
498              simError();
# Line 505 | Line 500 | void Integrator::constrainA(){
500  
501            rma = 1.0 / atoms[a]->getMass();
502            rmb = 1.0 / atoms[b]->getMass();
503 <          
503 >
504            gab = diffsq / ( 2.0 * ( rma + rmb ) * rpab );
505 +
506            dx = rxab * gab;
507            dy = ryab * gab;
508            dz = rzab * gab;
509  
510 <          pos[3*a+0] += rma * dx;
511 <          pos[3*a+1] += rma * dy;
512 <          pos[3*a+2] += rma * dz;
510 >          pos[ax] += rma * dx;
511 >          pos[ay] += rma * dy;
512 >          pos[az] += rma * dz;
513  
514 <          pos[3*b+0] -= rmb * dx;
515 <          pos[3*b+1] -= rmb * dy;
516 <          pos[3*b+2] -= rmb * dz;
514 >          pos[bx] -= rmb * dx;
515 >          pos[by] -= rmb * dy;
516 >          pos[bz] -= rmb * dz;
517  
518            dx = dx / dt;
519            dy = dy / dt;
520            dz = dz / dt;
521  
522 <          vel[3*a+0] += rma * dx;
523 <          vel[3*a+1] += rma * dy;
524 <          vel[3*a+2] += rma * dz;
522 >          vel[ax] += rma * dx;
523 >          vel[ay] += rma * dy;
524 >          vel[az] += rma * dz;
525  
526 <          vel[3*b+0] -= rmb * dx;
527 <          vel[3*b+1] -= rmb * dy;
528 <          vel[3*b+2] -= rmb * dz;
526 >          vel[bx] -= rmb * dx;
527 >          vel[by] -= rmb * dy;
528 >          vel[bz] -= rmb * dz;
529  
530            moving[a] = 1;
531            moving[b] = 1;
# Line 564 | Line 560 | void Integrator::constrainB( void ){
560    int done;
561    double vxab, vyab, vzab;
562    double rxab, ryab, rzab;
563 <  int a, b;
563 >  int a, b, ax, ay, az, bx, by, bz;
564    double rma, rmb;
565    double dx, dy, dz;
566    double rabsq, pabsq, rvab;
# Line 581 | Line 577 | void Integrator::constrainB( void ){
577    iteration = 0;
578    while( !done && (iteration < maxIteration ) ){
579  
580 +    done = 1;
581 +
582      for(i=0; i<nConstrained; i++){
583        
584        a = constrainedA[i];
585        b = constrainedB[i];
586  
587 +      ax = (a*3) + 0;
588 +      ay = (a*3) + 1;
589 +      az = (a*3) + 2;
590 +
591 +      bx = (b*3) + 0;
592 +      by = (b*3) + 1;
593 +      bz = (b*3) + 2;
594 +
595        if( moved[a] || moved[b] ){
596          
597 <        vxab = vel[3*a+0] - vel[3*b+0];
598 <        vyab = vel[3*a+1] - vel[3*b+1];
599 <        vzab = vel[3*a+2] - vel[3*b+2];
597 >        vxab = vel[ax] - vel[bx];
598 >        vyab = vel[ay] - vel[by];
599 >        vzab = vel[az] - vel[bz];
600  
601 <        rxab = pos[3*a+0] - pos[3*b+0];
602 <        ryab = pos[3*a+1] - pos[3*b+1];
603 <        rzab = pos[3*a+2] - pos[3*b+2];
601 >        rxab = pos[ax] - pos[bx];
602 >        ryab = pos[ay] - pos[by];
603 >        rzab = pos[az] - pos[bz];
604          
605 +
606          rxab = rxab - info->box_x * copysign(1, rxab)
607 <          * int( fabs(rxab) / info->box_x + 0.5);
607 >          * (int)( fabs(rxab / info->box_x) + 0.5);
608          ryab = ryab - info->box_y * copysign(1, ryab)
609 <          * int( fabs(ryab) / info->box_y + 0.5);
609 >          * (int)( fabs(ryab / info->box_y) + 0.5);
610          rzab = rzab - info->box_z * copysign(1, rzab)
611 <          * int( fabs(rzab) / info->box_z + 0.5);
612 <
611 >          * (int)( fabs(rzab / info->box_z) + 0.5);
612 >        
613          rma = 1.0 / atoms[a]->getMass();
614          rmb = 1.0 / atoms[b]->getMass();
615  
# Line 616 | Line 623 | void Integrator::constrainB( void ){
623            dy = ryab * gab;
624            dz = rzab * gab;
625            
626 <          vel[3*a+0] += rma * dx;
627 <          vel[3*a+1] += rma * dy;
628 <          vel[3*a+2] += rma * dz;
626 >          vel[ax] += rma * dx;
627 >          vel[ay] += rma * dy;
628 >          vel[az] += rma * dz;
629  
630 <          vel[3*b+0] -= rmb * dx;
631 <          vel[3*b+1] -= rmb * dy;
632 <          vel[3*b+2] -= rmb * dz;
630 >          vel[bx] -= rmb * dx;
631 >          vel[by] -= rmb * dy;
632 >          vel[bz] -= rmb * dz;
633            
634            moving[a] = 1;
635            moving[b] = 1;

Diff Legend

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