| 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(); | 
| 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 | 
  | 
  | 
| 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; |