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 556 by mmeineke, Wed Jun 18 22:20:12 2003 UTC

# Line 10 | Line 10
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 155 | Line 174 | void Symplectic::integrate( void ){
174    double A[3][3]; // the rotation matrix
175    double press[9];
176  
158 –  int time;
159 –
177    double dt          = info->dt;
178    double runTime     = info->run_time;
179    double sampleTime  = info->sampleTime;
180    double statusTime  = info->statusTime;
181    double thermalTime = info->thermalTime;
182  
183 <  int n_loops  = (int)( runTime / dt );
184 <  int sample_n = (int)( sampleTime / dt );
185 <  int status_n = (int)( statusTime / dt );
186 <  int vel_n    = (int)( thermalTime / dt );
183 >  double currSample;
184 >  double currThermal;
185 >  double currStatus;
186 >  double currTime;
187  
188    int calcPot, calcStress;
189    int isError;
# Line 191 | Line 208 | void Symplectic::integrate( void ){
208    dump_out->writeDump( 0.0 );
209    e_out->writeStat( 0.0 );
210    
211 <  calcPot = 0;
211 >  calcPot     = 0;
212 >  calcStress  = 0;
213 >  currSample  = sampleTime;
214 >  currThermal = thermalTime;
215 >  currStatus  = statusTime;
216 >  currTime    = 0.0;;
217  
218 <  for( tl=0; tl<nLoops; tl++){
218 >  while( currTime < runTime ){
219  
220 +    if( (currTime+dt) >= currStatus ){
221 +      calcPot = 1;
222 +      calcStress = 1;
223 +    }
224 +    
225      integrateStep( calcPot, calcStress );
226        
227 <    time = tl + 1;
228 <    
227 >    currTime += dt;
228 >
229      if( info->setTemp ){
230 <      if( !(time % vel_n) ) tStats->velocitize();
230 >      if( currTime >= currThermal ){
231 >        tStats->velocitize();
232 >        currThermal += thermalTime;
233 >      }
234      }
235 <    if( !(time % sample_n) ) dump_out->writeDump( time * dt );
236 <    if( !((time+1) % status_n) ) {
237 <      calcPot = 1;
238 <      calcStress = 1;
235 >
236 >    if( currTime >= currSample ){
237 >      dump_out->writeDump( currTime );
238 >      currSample += sampleTime;
239      }
240 <    if( !(time % status_n) ){
240 >
241 >    if( currTime >= currStatus ){
242        e_out->writeStat( time * dt );
243        calcPot = 0;
244 <      if (!strcasecmp(info->ensemble, "NPT")) calcStress = 1;
245 <      else calcStress = 0;
246 <    }    
216 <
217 <  
244 >      calcStress = 0;
245 >      currStatus += statusTime;
246 >    }
247    }
248  
249    dump_out->writeFinal();
# Line 223 | Line 252 | void Symplectic::integrate( void ){
252    delete e_out;
253   }
254  
255 + void Integrator::integrateStep( int calcPot, int calcStress ){
256  
257 < void Symplectic::moveA( void ){
257 >  // Position full step, and velocity half step
258 >
259 >  preMove();
260 >  moveA();
261 >  if( nConstrained ) constrainA();
262 >
263 >  // calc forces
264 >
265 >  myFF->doForces(calcPot,calcStress);
266 >
267 >  // finish the velocity  half step
268    
269 +  moveB();
270 +  if( nConstrained ) constrainB();
271 +  
272 + }
273 +
274 +
275 + void Integrator::moveA( void ){
276 +  
277    int i,j,k;
278    int atomIndex, aMatIndex;
279    DirectionalAtom* dAtom;
# Line 340 | Line 388 | void Integrator::moveB( void ){
388  
389   }
390  
391 + void Integrator::preMove( void ){
392 +  int i;
393  
394 +  if( nConstrained ){
395 +    if( oldAtoms != nAtoms ){
396 +      
397 +      // save oldAtoms to check for lode balanceing later on.
398 +      
399 +      oldAtoms = nAtoms;
400 +      
401 +      delete[] moving;
402 +      delete[] moved;
403 +      delete[] oldPos;
404 +      
405 +      moving = new int[nAtoms];
406 +      moved  = new int[nAtoms];
407 +      
408 +      oldPos = new double[nAtoms*3];
409 +    }
410 +  
411 +    for(i=0; i<(nAtoms*3); i++) oldPos[i] = pos[i];
412 +  }
413 + }  
414 +
415   void Integrator::constrainA(){
416  
417 +  int i,j,k;
418 +  int done;
419 +  double pxab, pyab, pzab;
420 +  double rxab, ryab, rzab;
421 +  int a, b;
422 +  double rma, rmb;
423 +  double dx, dy, dz;
424 +  double rabsq, pabsq, rpabsq;
425 +  double diffsq;
426 +  double gab;
427 +  int iteration;
428 +
429 +
430    
431 +  for( i=0; i<nAtoms; i++){
432 +    
433 +    moving[i] = 0;
434 +    moved[i]  = 1;
435 +  }
436 +  
437 +  
438 +  iteration = 0;
439 +  done = 0;
440 +  while( !done && (iteration < maxIteration )){
441  
442 +    done = 1;
443 +    for(i=0; i<nConstrained; i++){
444  
445 +      a = constrainedA[i];
446 +      b = constrainedB[i];
447 +    
448 +      if( moved[a] || moved[b] ){
449 +        
450 +        pxab = pos[3*a+0] - pos[3*b+0];
451 +        pyab = pos[3*a+1] - pos[3*b+1];
452 +        pzab = pos[3*a+2] - pos[3*b+2];
453 +
454 +        //periodic boundary condition
455 +        pxab = pxab - info->box_x * copysign(1, pxab)
456 +          * int(pxab / info->box_x + 0.5);
457 +        pyab = pyab - info->box_y * copysign(1, pyab)
458 +          * int(pyab / info->box_y + 0.5);
459 +        pzab = pzab - info->box_z * copysign(1, pzab)
460 +          * int(pzab / info->box_z + 0.5);
461 +      
462 +        pabsq = pxab * pxab + pyab * pyab + pzab * pzab;
463 +        rabsq = constraintedDsqr[i];
464 +        diffsq = pabsq - rabsq;
465 +
466 +        // the original rattle code from alan tidesley
467 +        if (fabs(diffsq) > tol*rabsq*2) {
468 +          rxab = oldPos[3*a+0] - oldPos[3*b+0];
469 +          ryab = oldPos[3*a+1] - oldPos[3*b+1];
470 +          rzab = oldPos[3*a+2] - oldPos[3*b+2];
471 +
472 +          rxab = rxab - info->box_x * copysign(1, rxab)
473 +            * int(rxab / info->box_x + 0.5);
474 +          ryab = ryab - info->box_y * copysign(1, ryab)
475 +            * int(ryab / info->box_y + 0.5);
476 +          rzab = rzab - info->box_z * copysign(1, rzab)
477 +            * int(rzab / info->box_z + 0.5);
478 +
479 +          rpab = rxab * pxab + ryab * pyab + rzab * pzab;
480 +          rpabsq = rpab * rpab;
481 +
482 +
483 +          if (rpabsq < (rabsq * -diffsq)){
484 + #ifdef IS_MPI
485 +            a = atoms[a]->getGlobalIndex();
486 +            b = atoms[b]->getGlobalIndex();
487 + #endif //is_mpi
488 +            sprintf( painCave.errMsg,
489 +                     "Constraint failure in constrainA at atom %d and %d\n.",
490 +                     a, b );
491 +            painCave.isFatal = 1;
492 +            simError();
493 +          }
494 +
495 +          rma = 1.0 / atoms[a]->getMass();
496 +          rmb = 1.0 / atoms[b]->getMass();
497 +          
498 +          gab = diffsq / ( 2.0 * ( rma + rmb ) * rpab );
499 +          dx = rxab * gab;
500 +          dy = ryab * gab;
501 +          dz = rzab * gab;
502 +
503 +          pos[3*a+0] += rma * dx;
504 +          pos[3*a+1] += rma * dy;
505 +          pos[3*a+2] += rma * dz;
506 +
507 +          pos[3*b+0] -= rmb * dx;
508 +          pos[3*b+1] -= rmb * dy;
509 +          pos[3*b+2] -= rmb * dz;
510 +
511 +          dx = dx / dt;
512 +          dy = dy / dt;
513 +          dz = dz / dt;
514 +
515 +          vel[3*a+0] += rma * dx;
516 +          vel[3*a+1] += rma * dy;
517 +          vel[3*a+2] += rma * dz;
518 +
519 +          vel[3*b+0] -= rmb * dx;
520 +          vel[3*b+1] -= rmb * dy;
521 +          vel[3*b+2] -= rmb * dz;
522 +
523 +          moving[a] = 1;
524 +          moving[b] = 1;
525 +          done = 0;
526 +        }
527 +      }
528 +    }
529 +    
530 +    for(i=0; i<nAtoms; i++){
531 +      
532 +      moved[i] = moving[i];
533 +      moving[i] = 0;
534 +    }
535 +
536 +    iteration++;
537 +  }
538 +
539 +  if( !done ){
540 +
541 +    sprintf( painCae.errMsg,
542 +             "Constraint failure in constrainA, too many iterations: %d\n",
543 +             iterations );
544 +    painCave.isFatal = 1;
545 +    simError();
546 +  }
547 +
548   }
549  
550 + void Integrator::constrainB( void ){
551 +  
552 +  int i,j,k;
553 +  int done;
554 +  double vxab, vyab, vzab;
555 +  double rxab, ryab, rzab;
556 +  int a, b;
557 +  double rma, rmb;
558 +  double dx, dy, dz;
559 +  double rabsq, pabsq, rvab;
560 +  double diffsq;
561 +  double gab;
562 +  int iteration;
563  
564 +  for(i=0; i<nAtom; i++){
565 +    moving[i] = 0;
566 +    moved[i] = 1;
567 +  }
568  
569 +  done = 0;
570 +  while( !done && (iteration < maxIteration ) ){
571  
572 +    for(i=0; i<nConstrained; i++){
573 +      
574 +      a = constrainedA[i];
575 +      b = constrainedB[i];
576  
577 +      if( moved[a] || moved[b] ){
578 +        
579 +        vxab = vel[3*a+0] - vel[3*b+0];
580 +        vyab = vel[3*a+1] - vel[3*b+1];
581 +        vzab = vel[3*a+2] - vel[3*b+2];
582  
583 +        rxab = pos[3*a+0] - pos[3*b+0];q
584 +        ryab = pos[3*a+1] - pos[3*b+1];
585 +        rzab = pos[3*a+2] - pos[3*b+2];
586 +        
587 +        rxab = rxab - info->box_x * copysign(1, rxab)
588 +          * int(rxab / info->box_x + 0.5);
589 +        ryab = ryab - info->box_y * copysign(1, ryab)
590 +          * int(ryab / info->box_y + 0.5);
591 +        rzab = rzab - info->box_z * copysign(1, rzab)
592 +          * int(rzab / info->box_z + 0.5);
593  
594 +        rma = 1.0 / atoms[a]->getMass();
595 +        rmb = 1.0 / atoms[b]->getMass();
596  
597 +        rvab = rxab * vxab + ryab * vyab + rzab * vzab;
598 +          
599 +        gab = -rvab / ( ( rma + rmb ) * constraintsDsqr[i] );
600  
601 < void Symplectic::rotate( int axes1, int axes2, double angle, double ji[3],
601 >        if (fabs(gab) > tol) {
602 >          
603 >          dx = rxab * gab;
604 >          dy = ryab * gab;
605 >          dz = rzab * gab;
606 >          
607 >          vel[3*a+0] += rma * dx;
608 >          vel[3*a+1] += rma * dy;
609 >          vel[3*a+2] += rma * dz;
610 >
611 >          vel[3*b+0] -= rmb * dx;
612 >          vel[3*b+1] -= rmb * dy;
613 >          vel[3*b+2] -= rmb * dz;
614 >          
615 >          moving[a] = 1;
616 >          moving[b] = 1;
617 >          done = 0;
618 >        }
619 >      }
620 >    }
621 >
622 >    for(i=0; i<nAtoms; i++){
623 >      moved[i] = moving[i];
624 >      moving[i] = 0;
625 >    }
626 >    
627 >    iteration++;
628 >  }
629 >
630 >  if( !done ){
631 >
632 >  
633 >    sprintf( painCae.errMsg,
634 >             "Constraint failure in constrainB, too many iterations: %d\n",
635 >             iterations );
636 >    painCave.isFatal = 1;
637 >    simError();
638 >  }
639 >
640 > }
641 >
642 >
643 >
644 >
645 >
646 >
647 >
648 > void Integrator::rotate( int axes1, int axes2, double angle, double ji[3],
649                           double A[3][3] ){
650  
651    int i,j,k;

Diff Legend

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