# | Line 66 | Line 66 | void SimSetup::createSim( void ){ | |
---|---|---|
66 | ||
67 | MakeStamps *the_stamps; | |
68 | Globals* the_globals; | |
69 | + | ExtendedSystem* the_extendedsystem; |
70 | int i, j; | |
71 | ||
72 | // get the stamps and globals; | |
# | Line 80 | Line 81 | void SimSetup::createSim( void ){ | |
81 | // get the ones we know are there, yet still may need some work. | |
82 | n_components = the_globals->getNComponents(); | |
83 | strcpy( force_field, the_globals->getForceField() ); | |
84 | + | |
85 | + | // get the ensemble and set up an extended system if we need it: |
86 | strcpy( ensemble, the_globals->getEnsemble() ); | |
87 | + | if( !strcasecmp( ensemble, "NPT" ) ) { |
88 | + | the_extendedsystem = new ExtendedSystem( simnfo ); |
89 | + | the_extendedsystem->setTargetTemp(the_globals->getTargetTemp()); |
90 | + | if (the_globals->haveTargetPressure()) |
91 | + | the_extendedsystem->setTargetPressure(the_globals->getTargetPressure()); |
92 | + | else { |
93 | + | sprintf( painCave.errMsg, |
94 | + | "SimSetup error: If you use the constant pressure\n" |
95 | + | " ensemble, you must set targetPressure.\n" |
96 | + | " This was found in the BASS file.\n"); |
97 | + | painCave.isFatal = 1; |
98 | + | simError(); |
99 | + | } |
100 | + | |
101 | + | if (the_globals->haveTauThermostat()) |
102 | + | the_extendedsystem->setTauThermostat(the_globals->getTauThermostat()); |
103 | + | else if (the_globals->haveQmass()) |
104 | + | the_extendedsystem->setQmass(the_globals->getQmass()); |
105 | + | else { |
106 | + | sprintf( painCave.errMsg, |
107 | + | "SimSetup error: If you use one of the constant temperature\n" |
108 | + | " ensembles, you must set either tauThermostat or qMass.\n" |
109 | + | " Neither of these was found in the BASS file.\n"); |
110 | + | painCave.isFatal = 1; |
111 | + | simError(); |
112 | + | } |
113 | + | |
114 | + | if (the_globals->haveTauBarostat()) |
115 | + | the_extendedsystem->setTauBarostat(the_globals->getTauBarostat()); |
116 | + | else { |
117 | + | sprintf( painCave.errMsg, |
118 | + | "SimSetup error: If you use the constant pressure\n" |
119 | + | " ensemble, you must set tauBarostat.\n" |
120 | + | " This was found in the BASS file.\n"); |
121 | + | painCave.isFatal = 1; |
122 | + | simError(); |
123 | + | } |
124 | + | |
125 | + | } else if ( !strcasecmp( ensemble, "NVT") ) { |
126 | + | the_extendedsystem = new ExtendedSystem( simnfo ); |
127 | + | the_extendedsystem->setTargetTemp(the_globals->getTargetTemp()); |
128 | + | |
129 | + | if (the_globals->haveTauThermostat()) |
130 | + | the_extendedsystem->setTauThermostat(the_globals->getTauThermostat()); |
131 | + | else if (the_globals->haveQmass()) |
132 | + | the_extendedsystem->setQmass(the_globals->getQmass()); |
133 | + | else { |
134 | + | sprintf( painCave.errMsg, |
135 | + | "SimSetup error: If you use one of the constant temperature\n" |
136 | + | " ensembles, you must set either tauThermostat or qMass.\n" |
137 | + | " Neither of these was found in the BASS file.\n"); |
138 | + | painCave.isFatal = 1; |
139 | + | simError(); |
140 | + | } |
141 | + | |
142 | + | } else if ( !strcasecmp( ensemble, "NVE") ) { |
143 | + | } else { |
144 | + | sprintf( painCave.errMsg, |
145 | + | "SimSetup Warning. Unrecognized Ensemble -> %s, " |
146 | + | "reverting to NVE for this simulation.\n", |
147 | + | ensemble ); |
148 | + | painCave.isFatal = 0; |
149 | + | simError(); |
150 | + | strcpy( ensemble, "NVE" ); |
151 | + | } |
152 | strcpy( simnfo->ensemble, ensemble ); | |
153 | ||
154 | strcpy( simnfo->mixingRule, the_globals->getMixingRule() ); | |
155 | simnfo->usePBC = the_globals->getPBC(); | |
156 | ||
157 | < | |
158 | < | |
159 | < | if( !strcmp( force_field, "TraPPE" ) ) the_ff = new TraPPEFF(); |
160 | < | else if( !strcmp( force_field, "DipoleTest" ) ) the_ff = new DipoleTestFF(); |
161 | < | else if( !strcmp( force_field, "TraPPE_Ex" ) ) the_ff = new TraPPE_ExFF(); |
162 | < | else if( !strcmp( force_field, "LJ" ) ) the_ff = new LJ_FF(); |
157 | > | int usesDipoles = 0; |
158 | > | if( !strcmp( force_field, "TraPPE_Ex" ) ){ |
159 | > | the_ff = new TraPPE_ExFF(); |
160 | > | usesDipoles = 1; |
161 | > | } |
162 | > | else if( !strcasecmp( force_field, "LJ" ) ) the_ff = new LJ_FF(); |
163 | else{ | |
164 | sprintf( painCave.errMsg, | |
165 | "SimSetup Error. Unrecognized force field -> %s\n", | |
# | Line 242 | Line 310 | void SimSetup::createSim( void ){ | |
310 | ||
311 | globalIndex = mpiSim->divideLabor(); | |
312 | ||
245 | – | |
246 | – | |
313 | // set up the local variables | |
314 | ||
315 | int localMol, allMol; | |
316 | int local_atoms, local_bonds, local_bends, local_torsions, local_SRI; | |
317 | + | |
318 | + | int* mol2proc = mpiSim->getMolToProcMap(); |
319 | + | int* molCompType = mpiSim->getMolComponentType(); |
320 | ||
321 | allMol = 0; | |
322 | localMol = 0; | |
# | Line 259 | Line 328 | void SimSetup::createSim( void ){ | |
328 | ||
329 | for( j=0; j<components_nmol[i]; j++ ){ | |
330 | ||
331 | < | if( mpiSim->getMyMolStart() <= allMol && |
263 | < | allMol <= mpiSim->getMyMolEnd() ){ |
331 | > | if( mol2proc[j] == worldRank ){ |
332 | ||
333 | local_atoms += comp_stamps[i]->getNAtoms(); | |
334 | local_bonds += comp_stamps[i]->getNBonds(); | |
# | Line 279 | Line 347 | void SimSetup::createSim( void ){ | |
347 | if( local_atoms != simnfo->n_atoms ){ | |
348 | sprintf( painCave.errMsg, | |
349 | "SimSetup error: mpiSim's localAtom (%d) and SimSetup's" | |
350 | < | " localAtom (%d) are note equal.\n", |
350 | > | " localAtom (%d) are not equal.\n", |
351 | simnfo->n_atoms, | |
352 | local_atoms ); | |
353 | painCave.isFatal = 1; | |
# | Line 304 | Line 372 | void SimSetup::createSim( void ){ | |
372 | Atom::createArrays(simnfo->n_atoms); | |
373 | the_atoms = new Atom*[simnfo->n_atoms]; | |
374 | the_molecules = new Molecule[simnfo->n_mol]; | |
375 | + | int molIndex; |
376 | ||
377 | + | // initialize the molecule's stampID's |
378 | ||
379 | + | #ifdef IS_MPI |
380 | + | |
381 | + | |
382 | + | molIndex = 0; |
383 | + | for(i=0; i<mpiSim->getTotNmol(); i++){ |
384 | + | |
385 | + | if(mol2proc[i] == worldRank ){ |
386 | + | the_molecules[molIndex].setStampID( molCompType[i] ); |
387 | + | the_molecules[molIndex].setMyIndex( molIndex ); |
388 | + | molIndex++; |
389 | + | } |
390 | + | } |
391 | + | |
392 | + | #else // is_mpi |
393 | + | |
394 | + | molIndex = 0; |
395 | + | for(i=0; i<n_components; i++){ |
396 | + | for(j=0; j<components_nmol[i]; j++ ){ |
397 | + | the_molecules[molIndex].setStampID( i ); |
398 | + | the_molecules[molIndex].setMyIndex( molIndex ); |
399 | + | molIndex++; |
400 | + | } |
401 | + | } |
402 | + | |
403 | + | |
404 | + | #endif // is_mpi |
405 | + | |
406 | + | |
407 | if( simnfo->n_SRI ){ | |
408 | + | |
409 | Exclude::createArray(simnfo->n_SRI); | |
410 | the_excludes = new Exclude*[simnfo->n_SRI]; | |
411 | + | for( int ex=0; ex<simnfo->n_SRI; ex++) the_excludes[ex] = new Exclude(ex); |
412 | simnfo->globalExcludes = new int; | |
413 | < | simnfo->n_exclude = tot_SRI; |
413 | > | simnfo->n_exclude = simnfo->n_SRI; |
414 | } | |
415 | else{ | |
416 | ||
# | Line 326 | Line 426 | void SimSetup::createSim( void ){ | |
426 | // set the arrays into the SimInfo object | |
427 | ||
428 | simnfo->atoms = the_atoms; | |
429 | < | simnfo->sr_interactions = the_sris; |
429 | > | simnfo->molecules = the_molecules; |
430 | simnfo->nGlobalExcludes = 0; | |
431 | simnfo->excludes = the_excludes; | |
432 | ||
# | Line 383 | Line 483 | void SimSetup::createSim( void ){ | |
483 | ||
484 | the_ff->setSimInfo( simnfo ); | |
485 | ||
486 | < | makeAtoms(); |
486 | > | makeMolecules(); |
487 | simnfo->identArray = new int[simnfo->n_atoms]; | |
488 | + | simnfo->molMembershipArray = new int[simnfo->n_atoms]; |
489 | for(i=0; i<simnfo->n_atoms; i++){ | |
490 | simnfo->identArray[i] = the_atoms[i]->getIdent(); | |
491 | + | } |
492 | + | for(i=0; i< simnfo->n_mol; i++) { |
493 | + | the_molecules[i].atomicRollCall(simnfo->molMembershipArray); |
494 | } | |
495 | ||
392 | – | if( tot_bonds ){ |
393 | – | makeBonds(); |
394 | – | } |
395 | – | |
396 | – | if( tot_bends ){ |
397 | – | makeBends(); |
398 | – | } |
399 | – | |
400 | – | if( tot_torsions ){ |
401 | – | makeTorsions(); |
402 | – | } |
403 | – | |
404 | – | |
496 | if (the_globals->getUseRF() ) { | |
497 | simnfo->useReactionField = 1; | |
498 | ||
# | Line 443 | Line 534 | void SimSetup::createSim( void ){ | |
534 | } | |
535 | simnfo->dielectric = the_globals->getDielectric(); | |
536 | } else { | |
537 | < | if (simnfo->n_dipoles) { |
537 | > | if (usesDipoles) { |
538 | ||
539 | if( !the_globals->haveECR() ){ | |
540 | sprintf( painCave.errMsg, | |
541 | < | "SimSetup Warning: using default value of 1/2 the smallest" |
541 | > | "SimSetup Warning: using default value of 1/2 the smallest " |
542 | "box length for the electrostaticCutoffRadius.\n" | |
543 | "I hope you have a very fast processor!\n"); | |
544 | painCave.isFatal = 0; | |
# | Line 463 | Line 554 | void SimSetup::createSim( void ){ | |
554 | ||
555 | if( !the_globals->haveEST() ){ | |
556 | sprintf( painCave.errMsg, | |
557 | < | "SimSetup Warning: using default value of 5% of the" |
557 | > | "SimSetup Warning: using default value of 5%% of the " |
558 | "electrostaticCutoffRadius for the " | |
559 | "electrostaticSkinThickness\n" | |
560 | ); | |
# | Line 638 | Line 729 | void SimSetup::createSim( void ){ | |
729 | ||
730 | // new AllLong( simnfo ); | |
731 | ||
641 | – | if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo, the_ff ); |
642 | – | if( !strcmp( force_field, "DipoleTest" ) ) new Symplectic( simnfo, the_ff ); |
643 | – | if( !strcmp( force_field, "TraPPE_Ex" ) ) new Symplectic( simnfo, the_ff ); |
644 | – | if( !strcmp( force_field, "LJ" ) ) new Verlet( *simnfo, the_ff ); |
732 | ||
733 | < | |
733 | > | if( !strcmp( force_field, "TraPPE_Ex" ) ){ |
734 | > | new Symplectic(simnfo, the_ff, the_extendedsystem); |
735 | > | } |
736 | > | else if( !strcmp( force_field, "LJ" ) ){ |
737 | > | new Verlet( *simnfo, the_ff, the_extendedsystem ); |
738 | > | } |
739 | > | else { |
740 | > | std::cerr << "I'm a bug.\n"; |
741 | > | fprintf( stderr, "Ima bug. stderr %s\n", force_field); |
742 | > | } |
743 | > | #ifdef IS_MPI |
744 | > | mpiSim->mpiRefresh(); |
745 | > | #endif |
746 | ||
747 | // initialize the Fortran | |
748 | < | |
748 | > | |
749 | > | |
750 | simnfo->refreshSim(); | |
751 | ||
752 | if( !strcmp( simnfo->mixingRule, "standard") ){ | |
# | Line 683 | Line 783 | void SimSetup::makeMolecules( void ){ | |
783 | BondStamp* currentBond; | |
784 | BendStamp* currentBend; | |
785 | TorsionStamp* currentTorsion; | |
786 | + | |
787 | + | bond_pair* theBonds; |
788 | + | bend_set* theBends; |
789 | + | torsion_set* theTorsions; |
790 | + | |
791 | ||
792 | //init the forceField paramters | |
793 | ||
794 | the_ff->readParams(); | |
795 | ||
796 | ||
797 | < | // init the molecules |
797 | > | // init the atoms |
798 | ||
799 | + | double ux, uy, uz, u, uSqr; |
800 | + | |
801 | atomOffset = 0; | |
802 | excludeOffset = 0; | |
803 | for(i=0; i<simnfo->n_mol; i++){ | |
# | Line 707 | Line 814 | void SimSetup::makeMolecules( void ){ | |
814 | info.myExcludes = &the_excludes[excludeOffset]; | |
815 | info.myBonds = new Bond*[info.nBonds]; | |
816 | info.myBends = new Bend*[info.nBends]; | |
817 | < | info.myTorsions = new Torsions*[info.nTorsions]; |
817 | > | info.myTorsions = new Torsion*[info.nTorsions]; |
818 | ||
819 | theBonds = new bond_pair[info.nBonds]; | |
820 | theBends = new bend_set[info.nBends]; | |
# | Line 717 | Line 824 | void SimSetup::makeMolecules( void ){ | |
824 | ||
825 | for(j=0; j<info.nAtoms; j++){ | |
826 | ||
827 | < | currentAtom = theComponents[stampID]->getAtom( j ); |
827 | > | currentAtom = comp_stamps[stampID]->getAtom( j ); |
828 | if( currentAtom->haveOrientation() ){ | |
829 | ||
830 | dAtom = new DirectionalAtom(j + atomOffset); | |
# | Line 758 | Line 865 | void SimSetup::makeMolecules( void ){ | |
865 | theBonds[j].a = currentBond->getA() + atomOffset; | |
866 | theBonds[j].b = currentBond->getB() + atomOffset; | |
867 | ||
868 | < | exI = theBonds[i].a; |
869 | < | exJ = theBonds[i].b; |
868 | > | exI = theBonds[j].a; |
869 | > | exJ = theBonds[j].b; |
870 | ||
871 | // exclude_I must always be the smaller of the pair | |
872 | if( exI > exJ ){ | |
# | Line 775 | Line 882 | void SimSetup::makeMolecules( void ){ | |
882 | ||
883 | the_excludes[j+excludeOffset]->setPair( exI, exJ ); | |
884 | #else // isn't MPI | |
885 | + | |
886 | the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) ); | |
887 | #endif //is_mpi | |
888 | } | |
# | Line 790 | Line 898 | void SimSetup::makeMolecules( void ){ | |
898 | ||
899 | if( currentBend->haveExtras() ){ | |
900 | ||
901 | < | extras = current_bend->getExtras(); |
901 | > | extras = currentBend->getExtras(); |
902 | current_extra = extras; | |
903 | ||
904 | while( current_extra != NULL ){ | |
# | Line 812 | Line 920 | void SimSetup::makeMolecules( void ){ | |
920 | ||
921 | default: | |
922 | sprintf( painCave.errMsg, | |
923 | < | "SimSetup Error: ghostVectorSource was neiter a " |
923 | > | "SimSetup Error: ghostVectorSource was neither a " |
924 | "double nor an int.\n" | |
925 | "-->Bend[%d] in %s\n", | |
926 | j, comp_stamps[stampID]->getID() ); | |
# | Line 906 | Line 1014 | void SimSetup::makeMolecules( void ){ | |
1014 | ||
1015 | ||
1016 | the_molecules[i].initialize( info ); | |
1017 | + | |
1018 | + | |
1019 | atomOffset += info.nAtoms; | |
1020 | + | delete[] theBonds; |
1021 | + | delete[] theBends; |
1022 | + | delete[] theTorsions; |
1023 | } | |
1024 | ||
1025 | + | #ifdef IS_MPI |
1026 | + | sprintf( checkPointMsg, "all molecules initialized succesfully" ); |
1027 | + | MPIcheckPoint(); |
1028 | + | #endif // is_mpi |
1029 | + | |
1030 | // clean up the forcefield | |
1031 | the_ff->calcRcut(); | |
1032 | the_ff->cleanMe(); | |
1033 | + | |
1034 | } | |
1035 | ||
1036 | void SimSetup::initFromBass( void ){ |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |