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

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 378 by mmeineke, Fri Mar 21 17:42:12 2003 UTC vs.
Revision 407 by mmeineke, Wed Mar 26 20:22:02 2003 UTC

# Line 333 | Line 333 | void SimSetup::createSim( void ){
333  
334    // get some of the tricky things that may still be in the globals
335  
336  if( simnfo->n_dipoles ){
337
338    if( !the_globals->haveRRF() ){
339      sprintf( painCave.errMsg,
340               "SimSetup Error, system has dipoles, but no rRF was set.\n");
341      painCave.isFatal = 1;
342      simError();
343    }
344    if( !the_globals->haveDielectric() ){
345      sprintf( painCave.errMsg,
346               "SimSetup Error, system has dipoles, but no"
347               " dielectric was set.\n" );
348      painCave.isFatal = 1;
349      simError();
350    }
351
352    simnfo->rRF        = the_globals->getRRF();
353    simnfo->dielectric = the_globals->getDielectric();
354  }
355
356 #ifdef IS_MPI
357  strcpy( checkPointMsg, "rRf and dielectric check out" );
358  MPIcheckPoint();
359 #endif // is_mpi
336    
337    if( the_globals->haveBox() ){
338      simnfo->box_x = the_globals->getBox();
# Line 426 | Line 402 | void SimSetup::createSim( void ){
402    }
403  
404  
405 <
405 >  if (the_globals->getUseRF() ) {
406 >    simnfo->useReactionField = 1;
407 >  
408 >    if( !the_globals->haveECR() ){
409 >      sprintf( painCave.errMsg,
410 >               "SimSetup Warning: using default value of 1/2 the smallest "
411 >               "box length for the electrostaticCutoffRadius.\n"
412 >               "I hope you have a very fast processor!\n");
413 >      painCave.isFatal = 0;
414 >      simError();
415 >      double smallest;
416 >      smallest = simnfo->box_x;
417 >      if (simnfo->box_y <= smallest) smallest = simnfo->box_y;
418 >      if (simnfo->box_z <= smallest) smallest = simnfo->box_z;
419 >      simnfo->ecr = 0.5 * smallest;
420 >    } else {
421 >      simnfo->ecr        = the_globals->getECR();
422 >    }
423  
424 +    if( !the_globals->haveEST() ){
425 +      sprintf( painCave.errMsg,
426 +               "SimSetup Warning: using default value of 0.05 * the "
427 +               "electrostaticCutoffRadius for the electrostaticSkinThickness\n"
428 +               );
429 +      painCave.isFatal = 0;
430 +      simError();
431 +      simnfo->est = 0.05 * simnfo->ecr;
432 +    } else {
433 +      simnfo->est        = the_globals->getEST();
434 +    }
435 +    
436 +    if(!the_globals->haveDielectric() ){
437 +      sprintf( painCave.errMsg,
438 +               "SimSetup Error: You are trying to use Reaction Field without"
439 +               "setting a dielectric constant!\n"
440 +               );
441 +      painCave.isFatal = 1;
442 +      simError();
443 +    }
444 +    simnfo->dielectric = the_globals->getDielectric();  
445 +  } else {
446 +    if (simnfo->n_dipoles) {
447 +      
448 +      if( !the_globals->haveECR() ){
449 +        sprintf( painCave.errMsg,
450 +                 "SimSetup Warning: using default value of 1/2 the smallest"
451 +                 "box length for the electrostaticCutoffRadius.\n"
452 +                 "I hope you have a very fast processor!\n");
453 +        painCave.isFatal = 0;
454 +        simError();
455 +        double smallest;
456 +        smallest = simnfo->box_x;
457 +        if (simnfo->box_y <= smallest) smallest = simnfo->box_y;
458 +        if (simnfo->box_z <= smallest) smallest = simnfo->box_z;
459 +        simnfo->ecr = 0.5 * smallest;
460 +      } else {
461 +        simnfo->ecr        = the_globals->getECR();
462 +      }
463 +      
464 +      if( !the_globals->haveEST() ){
465 +        sprintf( painCave.errMsg,
466 +                 "SimSetup Warning: using default value of 5% of the"
467 +                 "electrostaticCutoffRadius for the "
468 +                 "electrostaticSkinThickness\n"
469 +                 );
470 +        painCave.isFatal = 0;
471 +        simError();
472 +        simnfo->est = 0.05 * simnfo->ecr;
473 +      } else {
474 +        simnfo->est        = the_globals->getEST();
475 +      }
476 +    }
477 +  }  
478  
479 + #ifdef IS_MPI
480 +  strcpy( checkPointMsg, "electrostatic parameters check out" );
481 +  MPIcheckPoint();
482 + #endif // is_mpi
483  
484   if( the_globals->haveInitialConfig() ){
485  
# Line 620 | Line 671 | void SimSetup::makeAtoms( void ){
671   #endif // is_mpi
672   }
673  
674 +
675 + void SimSetup::makeMolecules( void ){
676 +
677 +  int i, j, exI, exJ, tempEx, stampID, atomOffset;
678 +  molInit info;
679 +  DirectionalAtom* dAtom;
680 +  AtomStamp* currentAtom;
681 +  BondStamp* currentBond;
682 +  BendStamp* currentBend;
683 +  TorsionStamp* currentTorsion;
684 +  
685 +  //init the forceField paramters
686 +
687 +  the_ff->readParams();
688 +
689 +  
690 +  // init the molecules
691 +
692 +  atomOffset = 0;
693 +  for(i=0; i<simnfo->n_mol; i++){
694 +    
695 +    stampID = the_molecules[i].getStampID();
696 +
697 +    info.nAtoms    = comp_stamps[stampID]->getNAtoms();
698 +    info.nBonds    = comp_stamps[stampID]->getNBonds();
699 +    info.nBends    = comp_stamps[stampID]->getNBends();
700 +    info.nTorsions = comp_stamps[stampID]->getNTorsions();
701 +    
702 +    info.myAtoms = &the_atoms[atomOffset];
703 +    info.myBonds = new Bond*[info.nBonds];
704 +    info.myBends = new Bend*[info.nBends];
705 +    info.myTorsions = new Torsions*[info.nTorsions];
706 +
707 +    theBonds = new bond_pair[info.nBonds];
708 +    theBends = new bend_set[info.nBends];
709 +    theTorsions = new torsion_set[info.nTorsions];
710 +    
711 +    // make the Atoms
712 +    
713 +    for(j=0; j<info.nAtoms; j++){
714 +      
715 +      currentAtom = theComponents[stampID]->getAtom( j );
716 +      if( currentAtom->haveOrientation() ){
717 +        
718 +        dAtom = new DirectionalAtom(j + atomOffset);
719 +        simnfo->n_oriented++;
720 +        info.myAtoms[j] = dAtom;
721 +        
722 +        ux = currentAtom->getOrntX();
723 +        uy = currentAtom->getOrntY();
724 +        uz = currentAtom->getOrntZ();
725 +        
726 +        uSqr = (ux * ux) + (uy * uy) + (uz * uz);
727 +        
728 +        u = sqrt( uSqr );
729 +        ux = ux / u;
730 +        uy = uy / u;
731 +        uz = uz / u;
732 +        
733 +        dAtom->setSUx( ux );
734 +        dAtom->setSUy( uy );
735 +        dAtom->setSUz( uz );
736 +      }
737 +      else{
738 +        info.myAtoms[j] = new GeneralAtom(j + atomOffset);
739 +      }
740 +      info.myAtoms[j]->setType( currentAtom->getType() );
741 +    
742 + #ifdef IS_MPI
743 +      
744 +      info.myAtoms[j]->setGlobalIndex( globalIndex[j+atomOffset] );
745 +      
746 + #endif // is_mpi
747 +    }
748 +    
749 +    // make the bonds
750 +    for(j=0; j<nBonds; j++){
751 +      
752 +      currentBond = comp_stamps[stampID]->getBond( j );
753 +      theBonds[j].a = currentBond->getA() + atomOffset;
754 +      theBonds[j].b = currentBond->getB() + atomOffset;
755 +
756 +      exI = theBonds[i].a;
757 +      exJ = theBonds[i].b;
758 +
759 +      // exclude_I must always be the smaller of the pair
760 +      if( exI > exJ ){
761 +        tempEx = exI;
762 +        exI = exJ;
763 +        exJ = tempEx;
764 +      }
765 + #ifdef IS_MPI
766 +      
767 +      the_excludes[index*2] =    
768 +        the_atoms[exI]->getGlobalIndex() + 1;
769 +      the_excludes[index*2 + 1] =
770 +        the_atoms[exJ]->getGlobalIndex() + 1;
771 +      
772 + #else  // isn't MPI
773 +      
774 +      the_excludes[index*2] =     exI + 1;
775 +      the_excludes[index*2 + 1] = exJ + 1;
776 +      // fortran index from 1 (hence the +1 in the indexing)
777 +
778 + #endif  //is_mpi
779 +    
780 +  }
781 +
782 +
783 +
784 +
785 +
786 +
787 +
788 +
789 +
790 +
791 +
792 +
793 +
794 +
795 +
796   void SimSetup::makeAtoms( void ){
797  
798    int i, j, k, index;

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines