# | 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::createSim( 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; |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |