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