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

Comparing trunk/OOPSE/libmdtools/Symplectic.cpp (file contents):
Revision 542 by mmeineke, Fri May 30 21:31:48 2003 UTC vs.
Revision 548 by mmeineke, Wed Jun 4 21:06:45 2003 UTC

# Line 10 | Line 10 | Symplectic::Symplectic( SimInfo* theInfo, ForceFields*
10   #include "simError.h"
11  
12  
13 < Symplectic::Symplectic( SimInfo* theInfo, ForceFields* the_ff ){
13 > Integrator::Integrator( SimInfo* theInfo, ForceFields* the_ff ){
14    
15    info = theInfo;
16    myFF = the_ff;
# Line 24 | Line 24 | Symplectic::Symplectic( SimInfo* theInfo, ForceFields*
24    if( info->the_integrator != NULL ) delete info->the_integrator;
25    info->the_integrator = this;
26  
27 +  nAtoms = info->n_atoms;
28 +
29    // check for constraints
30    
31 <  constrainedI = NULL;
32 <  constrainedJ = NULL;
31 >  constrainedA    = NULL;
32 >  constrainedB    = NULL;
33    constrainedDsqr = NULL;
34 +  moving          = NULL;
35 +  moved           = NULL;
36 +  prePos          = NULL;
37 +  
38    nConstrained = 0;
39  
40    checkConstraints();
41   }
42  
43 < Symplectic::~Symplectic() {
43 > Integrator::~Integrator() {
44    
45    if( nConstrained ){
46 <    delete[] constrainedI;
47 <    delete[] constrainedJ;
46 >    delete[] constrainedA;
47 >    delete[] constrainedB;
48      delete[] constrainedDsqr;
49 +    delete[] moving;
50 +    delete[] moved;
51 +    delete[] prePos;
52 + k
53    }
54    
55   }
56  
57 < void Symplectic::checkConstraints( void ){
57 > void Integrator::checkConstraints( void ){
58  
59  
60    isConstrained = 0;
# Line 114 | Line 124 | void Symplectic::checkConstraints( void ){
124      
125      isConstrained = 1;
126  
127 <    if(constrainedI != NULL )    delete[] constrainedI;
128 <    if(constrainedJ != NULL )    delete[] constrainedJ;
127 >    if(constrainedA != NULL )    delete[] constrainedA;
128 >    if(constrainedB != NULL )    delete[] constrainedB;
129      if(constrainedDsqr != NULL ) delete[] constrainedDsqr;
130  
131 <    constrainedI =    new int[nConstrained];
132 <    constrainedJ =    new int[nConstrained];
131 >    constrainedA =    new int[nConstrained];
132 >    constrainedB =    new int[nConstrained];
133      constrainedDsqr = new double[nConstrained];
134      
135      for( int i = 0; i < nConstrained; i++){
136        
137 <      constrainedI[i] = temp_con[i].get_a();
138 <      constrainedJ[i] = temp_con[i].get_b();
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 +    // save oldAtoms to check for lode balanceing later on.
144 +    
145 +    oldAtoms = nAtoms;
146 +    
147 +    moving = new int[nAtoms];
148 +    moved  = new int[nAtoms];
149 +
150 +    prePos = new double[nAtoms*3];
151    }
152    
153    delete[] temp_con;
154   }
155  
156  
157 < void Symplectic::integrate( void ){
157 > void Integrator::integrate( void ){
158  
159    int i, j;                         // loop counters
140 –  int nAtoms = info->n_atoms; // the number of atoms
160    double kE = 0.0;                  // the kinetic energy  
161    double rot_kE;
162    double trans_kE;
# Line 224 | Line 243 | void Symplectic::moveA( void ){
243   }
244  
245  
246 < void Symplectic::moveA( void ){
246 > void Integrator::moveA( void ){
247    
248    int i,j,k;
249    int atomIndex, aMatIndex;
# Line 340 | Line 359 | void Integrator::moveB( void ){
359  
360   }
361  
362 + void Integrator::preMove( void ){
363 +  int i;
364  
365 +  if( nConstrained ){
366 +    if( oldAtoms != nAtoms ){
367 +      
368 +      // save oldAtoms to check for lode balanceing later on.
369 +      
370 +      oldAtoms = nAtoms;
371 +      
372 +      delete[] moving;
373 +      delete[] moved;
374 +      delete[] oldPos;
375 +      
376 +      moving = new int[nAtoms];
377 +      moved  = new int[nAtoms];
378 +      
379 +      oldPos = new double[nAtoms*3];
380 +    }
381 +  
382 +    for(i=0; i<(nAtoms*3); i++) oldPos[i] = pos[i];
383 +  }
384 + }  
385 +
386   void Integrator::constrainA(){
387  
388 +  int i,j,k;
389 +  int done;
390 +  double pxab, pyab, pzab;
391 +  double rxab, ryab, rzab;
392 +  int a, b;
393 +  double rma, rmb;
394 +  double dx, dy, dz;
395 +  double rabsq, pabsq, rpabsq;
396 +  double diffsq;
397 +  double gab;
398 +  int iteration;
399 +
400 +
401    
402 +  for( i=0; i<nAtoms; i++){
403 +    
404 +    moving[i] = 0;
405 +    moved[i]  = 1;
406 +  }
407 +  
408 +  
409 +  iteration = 0;
410 +  done = 0;
411 +  while( !done && (iteration < maxIteration )){
412 +
413 +    done = 1;
414 +    for(i=0; i<nConstrained; i++){
415 +
416 +      a = constrainedA[i];
417 +      b = constrainedB[i];
418 +    
419 +      if( moved[a] || moved[b] ){
420 +        
421 +        pxab = pos[3*a+0] - pos[3*b+0];
422 +        pyab = pos[3*a+1] - pos[3*b+1];
423 +        pzab = pos[3*a+2] - pos[3*b+2];
424 +
425 +        //periodic boundary condition
426 +        pxab = pxab - info->box_x * copysign(1, pxab)
427 +          * int(pxab / info->box_x + 0.5);
428 +        pyab = pyab - info->box_y * copysign(1, pyab)
429 +          * int(pyab / info->box_y + 0.5);
430 +        pzab = pzab - info->box_z * copysign(1, pzab)
431 +          * int(pzab / info->box_z + 0.5);
432 +      
433 +        pabsq = pxab * pxab + pyab * pyab + pzab * pzab;
434 +        rabsq = constraintedDsqr[i];
435 +        diffsq = pabsq - rabsq;
436 +
437 +        // the original rattle code from alan tidesley
438 +        if (fabs(diffsq) > tol*rabsq*2) {
439 +          rxab = oldPos[3*a+0] - oldPos[3*b+0];
440 +          ryab = oldPos[3*a+1] - oldPos[3*b+1];
441 +          rzab = oldPos[3*a+2] - oldPos[3*b+2];
442 +
443 +          rxab = rxab - info->box_x * copysign(1, rxab)
444 +            * int(rxab / info->box_x + 0.5);
445 +          ryab = ryab - info->box_y * copysign(1, ryab)
446 +            * int(ryab / info->box_y + 0.5);
447 +          rzab = rzab - info->box_z * copysign(1, rzab)
448 +            * int(rzab / info->box_z + 0.5);
449 +
450 +          rpab = rxab * pxab + ryab * pyab + rzab * pzab;
451 +          rpabsq = rpab * rpab;
452 +
453 +
454 +          if (rpabsq < (rabsq * -diffsq)){
455 + #ifdef IS_MPI
456 +            a = atoms[a]->getGlobalIndex();
457 +            b = atoms[b]->getGlobalIndex();
458 + #endif //is_mpi
459 +            sprintf( painCave.errMsg,
460 +                     "Constraint failure in constrainA at atom %d and %d\n.",
461 +                     a, b );
462 +            painCave.isFatal = 1;
463 +            simError();
464 +          }
465  
466 +          rma = 1.0 / atoms[a]->getMass();
467 +          rmb = 1.0 / atoms[b]->getMass();
468 +          
469 +          gab = diffsq / ( 2.0 * ( rma + rmb ) * rpab );
470 +          dx = rxab * gab;
471 +          dy = ryab * gab;
472 +          dz = rzab * gab;
473  
474 +          pos[3*a+0] += rma * dx;
475 +          pos[3*a+1] += rma * dy;
476 +          pos[3*a+2] += rma * dz;
477 +
478 +          pos[3*b+0] -= rmb * dx;
479 +          pos[3*b+1] -= rmb * dy;
480 +          pos[3*b+2] -= rmb * dz;
481 +
482 +          dx = dx / dt;
483 +          dy = dy / dt;
484 +          dz = dz / dt;
485 +
486 +          vel[3*a+0] += rma * dx;
487 +          vel[3*a+1] += rma * dy;
488 +          vel[3*a+2] += rma * dz;
489 +
490 +          vel[3*b+0] -= rmb * dx;
491 +          vel[3*b+1] -= rmb * dy;
492 +          vel[3*b+2] -= rmb * dz;
493 +
494 +          moving[a] = 1;
495 +          moving[b] = 1;
496 +          done = 0;
497 +        }
498 +      }
499 +    }
500 +    
501 +    for(i=0; i<nAtoms; i++){
502 +      
503 +      moved[i] = moving[i];
504 +      moving[i] = 0;
505 +    }
506 +
507 +    iteration++;
508 +  }
509 +
510 +  if( !done ){
511 +
512 +    sprintf( painCae.errMsg,
513 +             "Constraint failure in constrainA, too many iterations: %d\n",
514 +             iterations );
515 +    painCave.isFatal = 1;
516 +    simError();
517 +  }
518 +
519   }
520  
521 + void Integrator::constrainB( void ){
522 +  
523 +  int i,j,k;
524 +  int done;
525 +  double vxab, vyab, vzab;
526 +  double rxab, ryab, rzab;
527 +  int a, b;
528 +  double rma, rmb;
529 +  double dx, dy, dz;
530 +  double rabsq, pabsq, rvab;
531 +  double diffsq;
532 +  double gab;
533 +  int iteration;
534  
535 +  for(i=0; i<nAtom; i++){
536 +    moving[i] = 0;
537 +    moved[i] = 1;
538 +  }
539  
540 +  done = 0;
541 +  while( !done && (iteration < maxIteration ) ){
542  
543 +    for(i=0; i<nConstrained; i++){
544 +      
545 +      a = constrainedA[i];
546 +      b = constrainedB[i];
547  
548 +      if( moved[a] || moved[b] ){
549 +        
550 +        vxab = vel[3*a+0] - vel[3*b+0];
551 +        vyab = vel[3*a+1] - vel[3*b+1];
552 +        vzab = vel[3*a+2] - vel[3*b+2];
553  
554 +        rxab = pos[3*a+0] - pos[3*b+0];q
555 +        ryab = pos[3*a+1] - pos[3*b+1];
556 +        rzab = pos[3*a+2] - pos[3*b+2];
557 +        
558 +        rxab = rxab - info->box_x * copysign(1, rxab)
559 +          * int(rxab / info->box_x + 0.5);
560 +        ryab = ryab - info->box_y * copysign(1, ryab)
561 +          * int(ryab / info->box_y + 0.5);
562 +        rzab = rzab - info->box_z * copysign(1, rzab)
563 +          * int(rzab / info->box_z + 0.5);
564  
565 +        rma = 1.0 / atoms[a]->getMass();
566 +        rmb = 1.0 / atoms[b]->getMass();
567  
568 +        rvab = rxab * vxab + ryab * vyab + rzab * vzab;
569 +          
570 +        gab = -rvab / ( ( rma + rmb ) * constraintsDsqr[i] );
571  
572 < void Symplectic::rotate( int axes1, int axes2, double angle, double ji[3],
572 >        if (fabs(gab) > tol) {
573 >          
574 >          dx = rxab * gab;
575 >          dy = ryab * gab;
576 >          dz = rzab * gab;
577 >          
578 >          vel[3*a+0] += rma * dx;
579 >          vel[3*a+1] += rma * dy;
580 >          vel[3*a+2] += rma * dz;
581 >
582 >          vel[3*b+0] -= rmb * dx;
583 >          vel[3*b+1] -= rmb * dy;
584 >          vel[3*b+2] -= rmb * dz;
585 >          
586 >          moving[a] = 1;
587 >          moving[b] = 1;
588 >          done = 0;
589 >        }
590 >      }
591 >    }
592 >
593 >    for(i=0; i<nAtoms; i++){
594 >      moved[i] = moving[i];
595 >      moving[i] = 0;
596 >    }
597 >    
598 >    iteration++;
599 >  }
600 >
601 >  if( !done ){
602 >
603 >  
604 >    sprintf( painCae.errMsg,
605 >             "Constraint failure in constrainB, too many iterations: %d\n",
606 >             iterations );
607 >    painCave.isFatal = 1;
608 >    simError();
609 >  }
610 >
611 > }
612 >
613 >
614 >
615 >
616 >
617 >
618 >
619 > void Integrator::rotate( int axes1, int axes2, double angle, double ji[3],
620                           double A[3][3] ){
621  
622    int i,j,k;

Diff Legend

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