# | Line 307 | Line 307 | void SimSetup::createSim( void ){ | |
---|---|---|
307 | ||
308 | ||
309 | if( simnfo->n_SRI ){ | |
310 | < | the_sris = new SRI*[simnfo->n_SRI]; |
311 | < | the_excludes = new int[2 * simnfo->n_SRI]; |
310 | > | Exclude::createArray(simnfo->n_SRI); |
311 | > | the_excludes = new Exclude*[simnfo->n_SRI]; |
312 | simnfo->globalExcludes = new int; | |
313 | simnfo->n_exclude = tot_SRI; | |
314 | } | |
315 | else{ | |
316 | ||
317 | < | the_excludes = new int[2]; |
318 | < | the_excludes[0] = 0; |
319 | < | the_excludes[1] = 0; |
317 | > | Exclude::createArray( 1 ); |
318 | > | the_excludes = new Exclude*; |
319 | > | the_excludes[0] = new Exclude(0); |
320 | > | the_excludes[0]->setPair( 0,0 ); |
321 | simnfo->globalExcludes = new int; | |
322 | simnfo->globalExcludes[0] = 0; | |
323 | < | |
323 | < | simnfo->n_exclude = 1; |
323 | > | simnfo->n_exclude = 0; |
324 | } | |
325 | ||
326 | // set the arrays into the SimInfo object | |
# | 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, excludeOffset; |
678 | + | molInit info; |
679 | + | DirectionalAtom* dAtom; |
680 | + | LinkedAssign* extras; |
681 | + | LinkedAssign* current_extra; |
682 | + | AtomStamp* currentAtom; |
683 | + | BondStamp* currentBond; |
684 | + | BendStamp* currentBend; |
685 | + | TorsionStamp* currentTorsion; |
686 | + | |
687 | + | //init the forceField paramters |
688 | + | |
689 | + | the_ff->readParams(); |
690 | + | |
691 | + | |
692 | + | // init the molecules |
693 | + | |
694 | + | atomOffset = 0; |
695 | + | excludeOffset = 0; |
696 | + | for(i=0; i<simnfo->n_mol; i++){ |
697 | + | |
698 | + | stampID = the_molecules[i].getStampID(); |
699 | + | |
700 | + | info.nAtoms = comp_stamps[stampID]->getNAtoms(); |
701 | + | info.nBonds = comp_stamps[stampID]->getNBonds(); |
702 | + | info.nBends = comp_stamps[stampID]->getNBends(); |
703 | + | info.nTorsions = comp_stamps[stampID]->getNTorsions(); |
704 | + | info.nExcludes = info.nBonds + info.nBends + info.nTorsions; |
705 | + | |
706 | + | info.myAtoms = &the_atoms[atomOffset]; |
707 | + | info.myExcludes = &the_excludes[excludeOffset]; |
708 | + | info.myBonds = new Bond*[info.nBonds]; |
709 | + | info.myBends = new Bend*[info.nBends]; |
710 | + | info.myTorsions = new Torsions*[info.nTorsions]; |
711 | + | |
712 | + | theBonds = new bond_pair[info.nBonds]; |
713 | + | theBends = new bend_set[info.nBends]; |
714 | + | theTorsions = new torsion_set[info.nTorsions]; |
715 | + | |
716 | + | // make the Atoms |
717 | + | |
718 | + | for(j=0; j<info.nAtoms; j++){ |
719 | + | |
720 | + | currentAtom = theComponents[stampID]->getAtom( j ); |
721 | + | if( currentAtom->haveOrientation() ){ |
722 | + | |
723 | + | dAtom = new DirectionalAtom(j + atomOffset); |
724 | + | simnfo->n_oriented++; |
725 | + | info.myAtoms[j] = dAtom; |
726 | + | |
727 | + | ux = currentAtom->getOrntX(); |
728 | + | uy = currentAtom->getOrntY(); |
729 | + | uz = currentAtom->getOrntZ(); |
730 | + | |
731 | + | uSqr = (ux * ux) + (uy * uy) + (uz * uz); |
732 | + | |
733 | + | u = sqrt( uSqr ); |
734 | + | ux = ux / u; |
735 | + | uy = uy / u; |
736 | + | uz = uz / u; |
737 | + | |
738 | + | dAtom->setSUx( ux ); |
739 | + | dAtom->setSUy( uy ); |
740 | + | dAtom->setSUz( uz ); |
741 | + | } |
742 | + | else{ |
743 | + | info.myAtoms[j] = new GeneralAtom(j + atomOffset); |
744 | + | } |
745 | + | info.myAtoms[j]->setType( currentAtom->getType() ); |
746 | + | |
747 | + | #ifdef IS_MPI |
748 | + | |
749 | + | info.myAtoms[j]->setGlobalIndex( globalIndex[j+atomOffset] ); |
750 | + | |
751 | + | #endif // is_mpi |
752 | + | } |
753 | + | |
754 | + | // make the bonds |
755 | + | for(j=0; j<info.nBonds; j++){ |
756 | + | |
757 | + | currentBond = comp_stamps[stampID]->getBond( j ); |
758 | + | theBonds[j].a = currentBond->getA() + atomOffset; |
759 | + | theBonds[j].b = currentBond->getB() + atomOffset; |
760 | + | |
761 | + | exI = theBonds[i].a; |
762 | + | exJ = theBonds[i].b; |
763 | + | |
764 | + | // exclude_I must always be the smaller of the pair |
765 | + | if( exI > exJ ){ |
766 | + | tempEx = exI; |
767 | + | exI = exJ; |
768 | + | exJ = tempEx; |
769 | + | } |
770 | + | #ifdef IS_MPI |
771 | + | tempEx = exI; |
772 | + | exI = the_atoms[tempEx]->getGlobalIndex() + 1; |
773 | + | tempEx = exJ; |
774 | + | exJ = the_atoms[tempEx]->getGlobalIndex() + 1; |
775 | + | |
776 | + | the_excludes[j+excludeOffset]->setPair( exI, exJ ); |
777 | + | #else // isn't MPI |
778 | + | the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) ); |
779 | + | #endif //is_mpi |
780 | + | } |
781 | + | excludeOffset += info.nBonds; |
782 | + | |
783 | + | //make the bends |
784 | + | for(j=0; j<info.nBends; j++){ |
785 | + | |
786 | + | currentBend = comp_stamps[stampID]->getBend( j ); |
787 | + | theBends[j].a = currentBend->getA() + atomOffset; |
788 | + | theBends[j].b = currentBend->getB() + atomOffset; |
789 | + | theBends[j].c = currentBend->getC() + atomOffset; |
790 | + | |
791 | + | if( currentBend->haveExtras() ){ |
792 | + | |
793 | + | extras = current_bend->getExtras(); |
794 | + | current_extra = extras; |
795 | + | |
796 | + | while( current_extra != NULL ){ |
797 | + | if( !strcmp( current_extra->getlhs(), "ghostVectorSource" )){ |
798 | + | |
799 | + | switch( current_extra->getType() ){ |
800 | + | |
801 | + | case 0: |
802 | + | theBends[j].ghost = |
803 | + | current_extra->getInt() + atomOffset; |
804 | + | theBends[j].isGhost = 1; |
805 | + | break; |
806 | + | |
807 | + | case 1: |
808 | + | theBends[j].ghost = |
809 | + | (int)current_extra->getDouble() + atomOffset; |
810 | + | theBends[j].isGhost = 1; |
811 | + | break; |
812 | + | |
813 | + | default: |
814 | + | sprintf( painCave.errMsg, |
815 | + | "SimSetup Error: ghostVectorSource was neiter a " |
816 | + | "double nor an int.\n" |
817 | + | "-->Bend[%d] in %s\n", |
818 | + | j, comp_stamps[stampID]->getID() ); |
819 | + | painCave.isFatal = 1; |
820 | + | simError(); |
821 | + | } |
822 | + | } |
823 | + | |
824 | + | else{ |
825 | + | |
826 | + | sprintf( painCave.errMsg, |
827 | + | "SimSetup Error: unhandled bend assignment:\n" |
828 | + | " -->%s in Bend[%d] in %s\n", |
829 | + | current_extra->getlhs(), |
830 | + | j, comp_stamps[stampID]->getID() ); |
831 | + | painCave.isFatal = 1; |
832 | + | simError(); |
833 | + | } |
834 | + | |
835 | + | current_extra = current_extra->getNext(); |
836 | + | } |
837 | + | } |
838 | + | |
839 | + | if( !theBends[j].isGhost ){ |
840 | + | |
841 | + | exI = theBends[j].a; |
842 | + | exJ = theBends[j].c; |
843 | + | } |
844 | + | else{ |
845 | + | |
846 | + | exI = theBends[j].a; |
847 | + | exJ = theBends[j].b; |
848 | + | } |
849 | + | |
850 | + | // exclude_I must always be the smaller of the pair |
851 | + | if( exI > exJ ){ |
852 | + | tempEx = exI; |
853 | + | exI = exJ; |
854 | + | exJ = tempEx; |
855 | + | } |
856 | + | #ifdef IS_MPI |
857 | + | tempEx = exI; |
858 | + | exI = the_atoms[tempEx]->getGlobalIndex() + 1; |
859 | + | tempEx = exJ; |
860 | + | exJ = the_atoms[tempEx]->getGlobalIndex() + 1; |
861 | + | |
862 | + | the_excludes[j+excludeOffset]->setPair( exI, exJ ); |
863 | + | #else // isn't MPI |
864 | + | the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) ); |
865 | + | #endif //is_mpi |
866 | + | } |
867 | + | excludeOffset += info.nBends; |
868 | + | |
869 | + | for(j=0; j<info.nTorsions; j++){ |
870 | + | |
871 | + | currentTorsion = comp_stamps[stampID]->getTorsion( j ); |
872 | + | theTorsions[j].a = currentTorsion->getA() + atomOffset; |
873 | + | theTorsions[j].b = currentTorsion->getB() + atomOffset; |
874 | + | theTorsions[j].c = currentTorsion->getC() + atomOffset; |
875 | + | theTorsions[j].d = currentTorsion->getD() + atomOffset; |
876 | + | |
877 | + | exI = theTorsions[j].a; |
878 | + | exJ = theTorsions[j].d; |
879 | + | |
880 | + | // exclude_I must always be the smaller of the pair |
881 | + | if( exI > exJ ){ |
882 | + | tempEx = exI; |
883 | + | exI = exJ; |
884 | + | exJ = tempEx; |
885 | + | } |
886 | + | #ifdef IS_MPI |
887 | + | tempEx = exI; |
888 | + | exI = the_atoms[tempEx]->getGlobalIndex() + 1; |
889 | + | tempEx = exJ; |
890 | + | exJ = the_atoms[tempEx]->getGlobalIndex() + 1; |
891 | + | |
892 | + | the_excludes[j+excludeOffset]->setPair( exI, exJ ); |
893 | + | #else // isn't MPI |
894 | + | the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) ); |
895 | + | #endif //is_mpi |
896 | + | } |
897 | + | excludeOffset += info.nTorsions; |
898 | + | |
899 | + | |
900 | + | |
901 | + | |
902 | + | |
903 | + | |
904 | + | |
905 | + | |
906 | void SimSetup::makeAtoms( void ){ | |
907 | ||
908 | int i, j, k, index; |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |