# | Line 66 | Line 66 | void SimSetup::createSim( void ){ | |
---|---|---|
66 | ||
67 | MakeStamps *the_stamps; | |
68 | Globals* the_globals; | |
69 | < | int i, j; |
69 | > | ExtendedSystem* the_extendedsystem; |
70 | > | int i, j, k, globalAtomIndex; |
71 | ||
72 | // get the stamps and globals; | |
73 | the_stamps = stamps; | |
# | 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_Ex" ) ) the_ff = new TraPPE_ExFF(); |
160 | < | 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 228 | Line 298 | void SimSetup::createSim( void ){ | |
298 | simnfo->n_torsions = tot_torsions; | |
299 | simnfo->n_SRI = tot_SRI; | |
300 | simnfo->n_mol = tot_nmol; | |
301 | + | simnfo->molMembershipArray = new int[tot_atoms]; |
302 | ||
232 | – | |
303 | #ifdef IS_MPI | |
304 | ||
305 | // divide the molecules among processors here. | |
306 | ||
307 | mpiSim = new mpiSimulation( simnfo ); | |
308 | ||
239 | – | |
240 | – | |
309 | globalIndex = mpiSim->divideLabor(); | |
310 | ||
243 | – | |
244 | – | |
311 | // set up the local variables | |
312 | ||
313 | int localMol, allMol; | |
# | Line 256 | Line 322 | void SimSetup::createSim( void ){ | |
322 | local_bonds = 0; | |
323 | local_bends = 0; | |
324 | local_torsions = 0; | |
325 | + | globalAtomIndex = 0; |
326 | + | |
327 | + | |
328 | for( i=0; i<n_components; i++ ){ | |
329 | ||
330 | for( j=0; j<components_nmol[i]; j++ ){ | |
331 | ||
332 | < | if( mol2proc[j] == worldRank ){ |
332 | > | if( mol2proc[allMol] == worldRank ){ |
333 | ||
334 | local_atoms += comp_stamps[i]->getNAtoms(); | |
335 | local_bonds += comp_stamps[i]->getNBonds(); | |
# | Line 268 | Line 337 | void SimSetup::createSim( void ){ | |
337 | local_torsions += comp_stamps[i]->getNTorsions(); | |
338 | localMol++; | |
339 | } | |
340 | < | allMol++; |
340 | > | for (k = 0; k < comp_stamps[i]->getNAtoms(); k++) { |
341 | > | simnfo->molMembershipArray[globalAtomIndex] = allMol; |
342 | > | globalAtomIndex++; |
343 | > | } |
344 | > | |
345 | > | allMol++; |
346 | } | |
347 | } | |
348 | local_SRI = local_bonds + local_bends + local_torsions; | |
349 | ||
350 | + | if (worldRank != 0) { |
351 | + | for (i =0; i < tot_atoms; i++){ |
352 | + | std::cerr << "i = " << i << " molMembershipArray[i] = " << simnfo->molMembershipArray[i] << "\n"; |
353 | + | } |
354 | + | } |
355 | ||
356 | simnfo->n_atoms = mpiSim->getMyNlocal(); | |
357 | ||
# | Line 316 | Line 395 | void SimSetup::createSim( void ){ | |
395 | ||
396 | if(mol2proc[i] == worldRank ){ | |
397 | the_molecules[molIndex].setStampID( molCompType[i] ); | |
398 | + | the_molecules[molIndex].setMyIndex( molIndex ); |
399 | + | the_molecules[molIndex].setGlobalIndex( i ); |
400 | molIndex++; | |
401 | } | |
402 | } | |
# | Line 323 | Line 404 | void SimSetup::createSim( void ){ | |
404 | #else // is_mpi | |
405 | ||
406 | molIndex = 0; | |
407 | + | globalAtomIndex = 0; |
408 | for(i=0; i<n_components; i++){ | |
409 | for(j=0; j<components_nmol[i]; j++ ){ | |
410 | the_molecules[molIndex].setStampID( i ); | |
411 | + | the_molecules[molIndex].setMyIndex( molIndex ); |
412 | + | the_molecules[molIndex].setGlobalIndex( molIndex ); |
413 | + | for (k = 0; k < comp_stamps[i]->getNAtoms(); k++) { |
414 | + | simnfo->molMembershipArray[globalAtomIndex] = molIndex; |
415 | + | globalAtomIndex++; |
416 | + | } |
417 | molIndex++; | |
418 | } | |
419 | } | |
# | Line 335 | Line 423 | void SimSetup::createSim( void ){ | |
423 | ||
424 | ||
425 | if( simnfo->n_SRI ){ | |
426 | + | |
427 | Exclude::createArray(simnfo->n_SRI); | |
428 | the_excludes = new Exclude*[simnfo->n_SRI]; | |
429 | + | for( int ex=0; ex<simnfo->n_SRI; ex++) the_excludes[ex] = new Exclude(ex); |
430 | simnfo->globalExcludes = new int; | |
431 | < | simnfo->n_exclude = tot_SRI; |
431 | > | simnfo->n_exclude = simnfo->n_SRI; |
432 | } | |
433 | else{ | |
434 | ||
# | Line 354 | Line 444 | void SimSetup::createSim( void ){ | |
444 | // set the arrays into the SimInfo object | |
445 | ||
446 | simnfo->atoms = the_atoms; | |
447 | + | simnfo->molecules = the_molecules; |
448 | simnfo->nGlobalExcludes = 0; | |
449 | simnfo->excludes = the_excludes; | |
450 | ||
# | Line 457 | Line 548 | void SimSetup::createSim( void ){ | |
548 | } | |
549 | simnfo->dielectric = the_globals->getDielectric(); | |
550 | } else { | |
551 | < | if (simnfo->n_dipoles) { |
551 | > | if (usesDipoles) { |
552 | ||
553 | if( !the_globals->haveECR() ){ | |
554 | sprintf( painCave.errMsg, | |
555 | < | "SimSetup Warning: using default value of 1/2 the smallest" |
555 | > | "SimSetup Warning: using default value of 1/2 the smallest " |
556 | "box length for the electrostaticCutoffRadius.\n" | |
557 | "I hope you have a very fast processor!\n"); | |
558 | painCave.isFatal = 0; | |
# | Line 477 | Line 568 | void SimSetup::createSim( void ){ | |
568 | ||
569 | if( !the_globals->haveEST() ){ | |
570 | sprintf( painCave.errMsg, | |
571 | < | "SimSetup Warning: using default value of 5% of the" |
571 | > | "SimSetup Warning: using default value of 5%% of the " |
572 | "electrostaticCutoffRadius for the " | |
573 | "electrostaticSkinThickness\n" | |
574 | ); | |
# | Line 652 | Line 743 | void SimSetup::createSim( void ){ | |
743 | ||
744 | // new AllLong( simnfo ); | |
745 | ||
655 | – | if( !strcmp( force_field, "TraPPE_Ex" ) ) new Symplectic( simnfo, the_ff ); |
656 | – | if( !strcmp( force_field, "LJ" ) ) new Verlet( *simnfo, the_ff ); |
746 | ||
747 | + | if( !strcmp( force_field, "TraPPE_Ex" ) ){ |
748 | + | new Symplectic(simnfo, the_ff, the_extendedsystem); |
749 | + | } |
750 | + | else if( !strcmp( force_field, "LJ" ) ){ |
751 | + | new Verlet( *simnfo, the_ff, the_extendedsystem ); |
752 | + | } |
753 | + | else { |
754 | + | std::cerr << "I'm a bug.\n"; |
755 | + | fprintf( stderr, "Ima bug. stderr %s\n", force_field); |
756 | + | } |
757 | + | #ifdef IS_MPI |
758 | + | mpiSim->mpiRefresh(); |
759 | + | #endif |
760 | ||
659 | – | |
761 | // initialize the Fortran | |
762 | < | |
762 | > | |
763 | > | |
764 | simnfo->refreshSim(); | |
765 | ||
766 | if( !strcmp( simnfo->mixingRule, "standard") ){ | |
# | Line 777 | Line 879 | void SimSetup::makeMolecules( void ){ | |
879 | theBonds[j].a = currentBond->getA() + atomOffset; | |
880 | theBonds[j].b = currentBond->getB() + atomOffset; | |
881 | ||
882 | < | exI = theBonds[i].a; |
883 | < | exJ = theBonds[i].b; |
882 | > | exI = theBonds[j].a; |
883 | > | exJ = theBonds[j].b; |
884 | ||
885 | // exclude_I must always be the smaller of the pair | |
886 | if( exI > exJ ){ | |
# | Line 794 | Line 896 | void SimSetup::makeMolecules( void ){ | |
896 | ||
897 | the_excludes[j+excludeOffset]->setPair( exI, exJ ); | |
898 | #else // isn't MPI | |
899 | + | |
900 | the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) ); | |
901 | #endif //is_mpi | |
902 | } | |
# | Line 831 | Line 934 | void SimSetup::makeMolecules( void ){ | |
934 | ||
935 | default: | |
936 | sprintf( painCave.errMsg, | |
937 | < | "SimSetup Error: ghostVectorSource was neiter a " |
937 | > | "SimSetup Error: ghostVectorSource was neither a " |
938 | "double nor an int.\n" | |
939 | "-->Bend[%d] in %s\n", | |
940 | j, comp_stamps[stampID]->getID() ); | |
# | Line 925 | Line 1028 | void SimSetup::makeMolecules( void ){ | |
1028 | ||
1029 | ||
1030 | the_molecules[i].initialize( info ); | |
1031 | + | |
1032 | + | |
1033 | atomOffset += info.nAtoms; | |
1034 | delete[] theBonds; | |
1035 | delete[] theBends; | |
1036 | delete[] theTorsions; | |
1037 | } | |
1038 | + | |
1039 | + | #ifdef IS_MPI |
1040 | + | sprintf( checkPointMsg, "all molecules initialized succesfully" ); |
1041 | + | MPIcheckPoint(); |
1042 | + | #endif // is_mpi |
1043 | ||
1044 | // clean up the forcefield | |
1045 | the_ff->calcRcut(); | |
1046 | the_ff->cleanMe(); | |
1047 | + | |
1048 | } | |
1049 | ||
1050 | void SimSetup::initFromBass( void ){ |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |