# | Line 9 | Line 9 | |
---|---|---|
9 | #include "parse_me.h" | |
10 | #include "Integrator.hpp" | |
11 | #include "simError.h" | |
12 | + | #include "RigidBody.hpp" |
13 | + | #include "OOPSEMinimizer.hpp" |
14 | + | //#include "ConstraintElement.hpp" |
15 | + | //#include "ConstraintPair.hpp" |
16 | ||
17 | #ifdef IS_MPI | |
18 | #include "mpiBASS.h" | |
# | Line 24 | Line 28 | |
28 | #define NPTxyz_ENS 4 | |
29 | ||
30 | ||
31 | < | #define FF_DUFF 0 |
32 | < | #define FF_LJ 1 |
33 | < | #define FF_EAM 2 |
31 | > | #define FF_DUFF 0 |
32 | > | #define FF_LJ 1 |
33 | > | #define FF_EAM 2 |
34 | > | #define FF_H2O 3 |
35 | ||
36 | using namespace std; | |
37 | ||
38 | + | /** |
39 | + | * Check whether dividend is divisble by divisor or not |
40 | + | */ |
41 | + | bool isDivisible(double dividend, double divisor){ |
42 | + | double tolerance = 0.000001; |
43 | + | double quotient; |
44 | + | double diff; |
45 | + | int intQuotient; |
46 | + | |
47 | + | quotient = dividend / divisor; |
48 | + | |
49 | + | if (quotient < 0) |
50 | + | quotient = -quotient; |
51 | + | |
52 | + | intQuotient = int (quotient + tolerance); |
53 | + | |
54 | + | diff = fabs(fabs(dividend) - intQuotient * fabs(divisor)); |
55 | + | |
56 | + | if (diff <= tolerance) |
57 | + | return true; |
58 | + | else |
59 | + | return false; |
60 | + | } |
61 | + | |
62 | SimSetup::SimSetup(){ | |
63 | ||
64 | initSuspend = false; | |
# | Line 103 | Line 132 | void SimSetup::createSim(void){ | |
132 | ||
133 | sysObjectsCreation(); | |
134 | ||
135 | + | // check on the post processing info |
136 | + | |
137 | + | finalInfoCheck(); |
138 | + | |
139 | // initialize the system coordinates | |
140 | ||
141 | if ( !initSuspend ){ | |
# | Line 112 | Line 145 | void SimSetup::createSim(void){ | |
145 | info[0].currentTime = 0.0; | |
146 | } | |
147 | ||
115 | – | // check on the post processing info |
116 | – | |
117 | – | finalInfoCheck(); |
118 | – | |
148 | // make the output filenames | |
149 | ||
150 | makeOutNames(); | |
151 | < | |
123 | < | // make the integrator |
124 | < | |
125 | < | makeIntegrator(); |
126 | < | |
151 | > | |
152 | #ifdef IS_MPI | |
153 | mpiSim->mpiRefresh(); | |
154 | #endif | |
# | Line 131 | Line 156 | void SimSetup::createSim(void){ | |
156 | // initialize the Fortran | |
157 | ||
158 | initFortran(); | |
159 | + | |
160 | + | if (globals->haveMinimizer()) |
161 | + | // make minimizer |
162 | + | makeMinimizer(); |
163 | + | else |
164 | + | // make the integrator |
165 | + | makeIntegrator(); |
166 | + | |
167 | } | |
168 | ||
169 | ||
170 | void SimSetup::makeMolecules(void){ | |
171 | < | int k; |
172 | < | int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset; |
171 | > | int i, j, k; |
172 | > | int exI, exJ, exK, exL, slI, slJ; |
173 | > | int tempI, tempJ, tempK, tempL; |
174 | > | int molI, globalID; |
175 | > | int stampID, atomOffset, rbOffset, groupOffset; |
176 | molInit molInfo; | |
177 | DirectionalAtom* dAtom; | |
178 | + | RigidBody* myRB; |
179 | + | StuntDouble* mySD; |
180 | LinkedAssign* extras; | |
181 | LinkedAssign* current_extra; | |
182 | AtomStamp* currentAtom; | |
183 | BondStamp* currentBond; | |
184 | BendStamp* currentBend; | |
185 | TorsionStamp* currentTorsion; | |
186 | + | RigidBodyStamp* currentRigidBody; |
187 | + | CutoffGroupStamp* currentCutoffGroup; |
188 | + | CutoffGroup* myCutoffGroup; |
189 | + | int nCutoffGroups;// number of cutoff group of a molecule defined in mdl file |
190 | + | set<int> cutoffAtomSet; //atoms belong to cutoffgroup defined at mdl file |
191 | ||
192 | bond_pair* theBonds; | |
193 | bend_set* theBends; | |
194 | torsion_set* theTorsions; | |
195 | ||
196 | + | set<int> skipList; |
197 | ||
198 | + | double phi, theta, psi; |
199 | + | char* molName; |
200 | + | char rbName[100]; |
201 | + | |
202 | + | //ConstraintPair* consPair; //constraint pair |
203 | + | //ConstraintElement* consElement1; //first element of constraint pair |
204 | + | //ConstraintElement* consElement2; //second element of constraint pair |
205 | + | //int whichRigidBody; |
206 | + | //int consAtomIndex; //index of constraint atom in rigid body's atom array |
207 | + | //vector<pair<int, int> > jointAtoms; |
208 | //init the forceField paramters | |
209 | ||
210 | the_ff->readParams(); | |
211 | ||
158 | – | |
212 | // init the atoms | |
213 | ||
214 | < | double ux, uy, uz, u, uSqr; |
214 | > | int nMembers, nNew, rb1, rb2; |
215 | ||
216 | for (k = 0; k < nInfo; k++){ | |
217 | the_ff->setSimInfo(&(info[k])); | |
218 | ||
219 | + | #ifdef IS_MPI |
220 | + | info[k].globalGroupMembership = new int[mpiSim->getNAtomsGlobal()]; |
221 | + | for (i = 0; i < mpiSim->getNAtomsGlobal(); i++) |
222 | + | info[k].globalGroupMembership[i] = 0; |
223 | + | #else |
224 | + | info[k].globalGroupMembership = new int[info[k].n_atoms]; |
225 | + | for (i = 0; i < info[k].n_atoms; i++) |
226 | + | info[k].globalGroupMembership[i] = 0; |
227 | + | #endif |
228 | + | |
229 | atomOffset = 0; | |
230 | < | excludeOffset = 0; |
230 | > | groupOffset = 0; |
231 | > | |
232 | for (i = 0; i < info[k].n_mol; i++){ | |
233 | stampID = info[k].molecules[i].getStampID(); | |
234 | + | molName = comp_stamps[stampID]->getID(); |
235 | ||
236 | molInfo.nAtoms = comp_stamps[stampID]->getNAtoms(); | |
237 | molInfo.nBonds = comp_stamps[stampID]->getNBonds(); | |
238 | molInfo.nBends = comp_stamps[stampID]->getNBends(); | |
239 | molInfo.nTorsions = comp_stamps[stampID]->getNTorsions(); | |
240 | < | molInfo.nExcludes = molInfo.nBonds + molInfo.nBends + molInfo.nTorsions; |
240 | > | molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies(); |
241 | ||
242 | + | nCutoffGroups = comp_stamps[stampID]->getNCutoffGroups(); |
243 | + | |
244 | molInfo.myAtoms = &(info[k].atoms[atomOffset]); | |
178 | – | molInfo.myExcludes = &(info[k].excludes[excludeOffset]); |
179 | – | molInfo.myBonds = new Bond * [molInfo.nBonds]; |
180 | – | molInfo.myBends = new Bend * [molInfo.nBends]; |
181 | – | molInfo.myTorsions = new Torsion * [molInfo.nTorsions]; |
245 | ||
246 | + | if (molInfo.nBonds > 0) |
247 | + | molInfo.myBonds = new Bond*[molInfo.nBonds]; |
248 | + | else |
249 | + | molInfo.myBonds = NULL; |
250 | + | |
251 | + | if (molInfo.nBends > 0) |
252 | + | molInfo.myBends = new Bend*[molInfo.nBends]; |
253 | + | else |
254 | + | molInfo.myBends = NULL; |
255 | + | |
256 | + | if (molInfo.nTorsions > 0) |
257 | + | molInfo.myTorsions = new Torsion *[molInfo.nTorsions]; |
258 | + | else |
259 | + | molInfo.myTorsions = NULL; |
260 | + | |
261 | theBonds = new bond_pair[molInfo.nBonds]; | |
262 | theBends = new bend_set[molInfo.nBends]; | |
263 | theTorsions = new torsion_set[molInfo.nTorsions]; | |
264 | < | |
264 | > | |
265 | // make the Atoms | |
266 | ||
267 | for (j = 0; j < molInfo.nAtoms; j++){ | |
268 | currentAtom = comp_stamps[stampID]->getAtom(j); | |
269 | + | |
270 | if (currentAtom->haveOrientation()){ | |
271 | dAtom = new DirectionalAtom((j + atomOffset), | |
272 | info[k].getConfiguration()); | |
273 | info[k].n_oriented++; | |
274 | molInfo.myAtoms[j] = dAtom; | |
275 | ||
276 | < | ux = currentAtom->getOrntX(); |
277 | < | uy = currentAtom->getOrntY(); |
278 | < | uz = currentAtom->getOrntZ(); |
276 | > | // Directional Atoms have standard unit vectors which are oriented |
277 | > | // in space using the three Euler angles. We assume the standard |
278 | > | // unit vector was originally along the z axis below. |
279 | ||
280 | < | uSqr = (ux * ux) + (uy * uy) + (uz * uz); |
280 | > | phi = currentAtom->getEulerPhi() * M_PI / 180.0; |
281 | > | theta = currentAtom->getEulerTheta() * M_PI / 180.0; |
282 | > | psi = currentAtom->getEulerPsi()* M_PI / 180.0; |
283 | ||
284 | < | u = sqrt(uSqr); |
285 | < | ux = ux / u; |
205 | < | uy = uy / u; |
206 | < | uz = uz / u; |
207 | < | |
208 | < | dAtom->setSUx(ux); |
209 | < | dAtom->setSUy(uy); |
210 | < | dAtom->setSUz(uz); |
284 | > | dAtom->setUnitFrameFromEuler(phi, theta, psi); |
285 | > | |
286 | } | |
287 | else{ | |
213 | – | molInfo.myAtoms[j] = new GeneralAtom((j + atomOffset), |
214 | – | info[k].getConfiguration()); |
215 | – | } |
216 | – | molInfo.myAtoms[j]->setType(currentAtom->getType()); |
288 | ||
289 | < | #ifdef IS_MPI |
289 | > | molInfo.myAtoms[j] = new Atom((j + atomOffset), info[k].getConfiguration()); |
290 | ||
291 | < | molInfo.myAtoms[j]->setGlobalIndex(globalIndex[j + atomOffset]); |
291 | > | } |
292 | ||
293 | + | molInfo.myAtoms[j]->setType(currentAtom->getType()); |
294 | + | #ifdef IS_MPI |
295 | + | molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]); |
296 | #endif // is_mpi | |
297 | } | |
298 | ||
# | Line 228 | Line 302 | void SimSetup::makeMolecules(void){ | |
302 | theBonds[j].a = currentBond->getA() + atomOffset; | |
303 | theBonds[j].b = currentBond->getB() + atomOffset; | |
304 | ||
305 | < | exI = theBonds[j].a; |
306 | < | exJ = theBonds[j].b; |
305 | > | tempI = theBonds[j].a; |
306 | > | tempJ = theBonds[j].b; |
307 | ||
234 | – | // exclude_I must always be the smaller of the pair |
235 | – | if (exI > exJ){ |
236 | – | tempEx = exI; |
237 | – | exI = exJ; |
238 | – | exJ = tempEx; |
239 | – | } |
308 | #ifdef IS_MPI | |
309 | < | tempEx = exI; |
310 | < | exI = info[k].atoms[tempEx]->getGlobalIndex() + 1; |
311 | < | tempEx = exJ; |
312 | < | exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1; |
309 | > | exI = info[k].atoms[tempI]->getGlobalIndex() + 1; |
310 | > | exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1; |
311 | > | #else |
312 | > | exI = tempI + 1; |
313 | > | exJ = tempJ + 1; |
314 | > | #endif |
315 | ||
316 | < | info[k].excludes[j + excludeOffset]->setPair(exI, exJ); |
247 | < | #else // isn't MPI |
248 | < | |
249 | < | info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1)); |
250 | < | #endif //is_mpi |
316 | > | info[k].excludes->addPair(exI, exJ); |
317 | } | |
252 | – | excludeOffset += molInfo.nBonds; |
318 | ||
319 | //make the bends | |
320 | for (j = 0; j < molInfo.nBends; j++){ | |
# | Line 299 | Line 364 | void SimSetup::makeMolecules(void){ | |
364 | } | |
365 | } | |
366 | ||
367 | < | if (!theBends[j].isGhost){ |
368 | < | exI = theBends[j].a; |
369 | < | exJ = theBends[j].c; |
370 | < | } |
371 | < | else{ |
307 | < | exI = theBends[j].a; |
308 | < | exJ = theBends[j].b; |
309 | < | } |
310 | < | |
311 | < | // exclude_I must always be the smaller of the pair |
312 | < | if (exI > exJ){ |
313 | < | tempEx = exI; |
314 | < | exI = exJ; |
315 | < | exJ = tempEx; |
316 | < | } |
367 | > | if (theBends[j].isGhost) { |
368 | > | |
369 | > | tempI = theBends[j].a; |
370 | > | tempJ = theBends[j].b; |
371 | > | |
372 | #ifdef IS_MPI | |
373 | < | tempEx = exI; |
374 | < | exI = info[k].atoms[tempEx]->getGlobalIndex() + 1; |
375 | < | tempEx = exJ; |
376 | < | exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1; |
373 | > | exI = info[k].atoms[tempI]->getGlobalIndex() + 1; |
374 | > | exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1; |
375 | > | #else |
376 | > | exI = tempI + 1; |
377 | > | exJ = tempJ + 1; |
378 | > | #endif |
379 | > | info[k].excludes->addPair(exI, exJ); |
380 | ||
381 | < | info[k].excludes[j + excludeOffset]->setPair(exI, exJ); |
382 | < | #else // isn't MPI |
383 | < | info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1)); |
384 | < | #endif //is_mpi |
381 | > | } else { |
382 | > | |
383 | > | tempI = theBends[j].a; |
384 | > | tempJ = theBends[j].b; |
385 | > | tempK = theBends[j].c; |
386 | > | |
387 | > | #ifdef IS_MPI |
388 | > | exI = info[k].atoms[tempI]->getGlobalIndex() + 1; |
389 | > | exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1; |
390 | > | exK = info[k].atoms[tempK]->getGlobalIndex() + 1; |
391 | > | #else |
392 | > | exI = tempI + 1; |
393 | > | exJ = tempJ + 1; |
394 | > | exK = tempK + 1; |
395 | > | #endif |
396 | > | |
397 | > | info[k].excludes->addPair(exI, exK); |
398 | > | info[k].excludes->addPair(exI, exJ); |
399 | > | info[k].excludes->addPair(exJ, exK); |
400 | > | } |
401 | } | |
328 | – | excludeOffset += molInfo.nBends; |
402 | ||
403 | for (j = 0; j < molInfo.nTorsions; j++){ | |
404 | currentTorsion = comp_stamps[stampID]->getTorsion(j); | |
# | Line 334 | Line 407 | void SimSetup::makeMolecules(void){ | |
407 | theTorsions[j].c = currentTorsion->getC() + atomOffset; | |
408 | theTorsions[j].d = currentTorsion->getD() + atomOffset; | |
409 | ||
410 | < | exI = theTorsions[j].a; |
411 | < | exJ = theTorsions[j].d; |
410 | > | tempI = theTorsions[j].a; |
411 | > | tempJ = theTorsions[j].b; |
412 | > | tempK = theTorsions[j].c; |
413 | > | tempL = theTorsions[j].d; |
414 | ||
340 | – | // exclude_I must always be the smaller of the pair |
341 | – | if (exI > exJ){ |
342 | – | tempEx = exI; |
343 | – | exI = exJ; |
344 | – | exJ = tempEx; |
345 | – | } |
415 | #ifdef IS_MPI | |
416 | < | tempEx = exI; |
417 | < | exI = info[k].atoms[tempEx]->getGlobalIndex() + 1; |
418 | < | tempEx = exJ; |
419 | < | exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1; |
416 | > | exI = info[k].atoms[tempI]->getGlobalIndex() + 1; |
417 | > | exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1; |
418 | > | exK = info[k].atoms[tempK]->getGlobalIndex() + 1; |
419 | > | exL = info[k].atoms[tempL]->getGlobalIndex() + 1; |
420 | > | #else |
421 | > | exI = tempI + 1; |
422 | > | exJ = tempJ + 1; |
423 | > | exK = tempK + 1; |
424 | > | exL = tempL + 1; |
425 | > | #endif |
426 | ||
427 | < | info[k].excludes[j + excludeOffset]->setPair(exI, exJ); |
428 | < | #else // isn't MPI |
429 | < | info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1)); |
430 | < | #endif //is_mpi |
431 | < | } |
432 | < | excludeOffset += molInfo.nTorsions; |
427 | > | info[k].excludes->addPair(exI, exJ); |
428 | > | info[k].excludes->addPair(exI, exK); |
429 | > | info[k].excludes->addPair(exI, exL); |
430 | > | info[k].excludes->addPair(exJ, exK); |
431 | > | info[k].excludes->addPair(exJ, exL); |
432 | > | info[k].excludes->addPair(exK, exL); |
433 | > | } |
434 | ||
435 | + | |
436 | + | molInfo.myRigidBodies.clear(); |
437 | + | |
438 | + | for (j = 0; j < molInfo.nRigidBodies; j++){ |
439 | ||
440 | < | // send the arrays off to the forceField for init. |
440 | > | currentRigidBody = comp_stamps[stampID]->getRigidBody(j); |
441 | > | nMembers = currentRigidBody->getNMembers(); |
442 | ||
443 | + | // Create the Rigid Body: |
444 | + | |
445 | + | myRB = new RigidBody(); |
446 | + | |
447 | + | sprintf(rbName,"%s_RB_%d", molName, j); |
448 | + | myRB->setType(rbName); |
449 | + | |
450 | + | for (rb1 = 0; rb1 < nMembers; rb1++) { |
451 | + | |
452 | + | // molI is atom numbering inside this molecule |
453 | + | molI = currentRigidBody->getMember(rb1); |
454 | + | |
455 | + | // tempI is atom numbering on local processor |
456 | + | tempI = molI + atomOffset; |
457 | + | |
458 | + | // currentAtom is the AtomStamp (which we need for |
459 | + | // rigid body reference positions) |
460 | + | currentAtom = comp_stamps[stampID]->getAtom(molI); |
461 | + | |
462 | + | // When we add to the rigid body, add the atom itself and |
463 | + | // the stamp info: |
464 | + | |
465 | + | myRB->addAtom(info[k].atoms[tempI], currentAtom); |
466 | + | |
467 | + | // Add this atom to the Skip List for the integrators |
468 | + | #ifdef IS_MPI |
469 | + | slI = info[k].atoms[tempI]->getGlobalIndex(); |
470 | + | #else |
471 | + | slI = tempI; |
472 | + | #endif |
473 | + | skipList.insert(slI); |
474 | + | |
475 | + | } |
476 | + | |
477 | + | for(rb1 = 0; rb1 < nMembers - 1; rb1++) { |
478 | + | for(rb2 = rb1+1; rb2 < nMembers; rb2++) { |
479 | + | |
480 | + | tempI = currentRigidBody->getMember(rb1); |
481 | + | tempJ = currentRigidBody->getMember(rb2); |
482 | + | |
483 | + | // Some explanation is required here. |
484 | + | // Fortran indexing starts at 1, while c indexing starts at 0 |
485 | + | // Also, in parallel computations, the GlobalIndex is |
486 | + | // used for the exclude list: |
487 | + | |
488 | + | #ifdef IS_MPI |
489 | + | exI = molInfo.myAtoms[tempI]->getGlobalIndex() + 1; |
490 | + | exJ = molInfo.myAtoms[tempJ]->getGlobalIndex() + 1; |
491 | + | #else |
492 | + | exI = molInfo.myAtoms[tempI]->getIndex() + 1; |
493 | + | exJ = molInfo.myAtoms[tempJ]->getIndex() + 1; |
494 | + | #endif |
495 | + | |
496 | + | info[k].excludes->addPair(exI, exJ); |
497 | + | |
498 | + | } |
499 | + | } |
500 | + | |
501 | + | molInfo.myRigidBodies.push_back(myRB); |
502 | + | info[k].rigidBodies.push_back(myRB); |
503 | + | } |
504 | + | |
505 | + | |
506 | + | //create cutoff group for molecule |
507 | + | |
508 | + | cutoffAtomSet.clear(); |
509 | + | molInfo.myCutoffGroups.clear(); |
510 | + | |
511 | + | for (j = 0; j < nCutoffGroups; j++){ |
512 | + | |
513 | + | currentCutoffGroup = comp_stamps[stampID]->getCutoffGroup(j); |
514 | + | nMembers = currentCutoffGroup->getNMembers(); |
515 | + | |
516 | + | myCutoffGroup = new CutoffGroup(); |
517 | + | |
518 | + | #ifdef IS_MPI |
519 | + | myCutoffGroup->setGlobalIndex(globalGroupIndex[groupOffset]); |
520 | + | #else |
521 | + | myCutoffGroup->setGlobalIndex(groupOffset); |
522 | + | #endif |
523 | + | |
524 | + | for (int cg = 0; cg < nMembers; cg++) { |
525 | + | |
526 | + | // molI is atom numbering inside this molecule |
527 | + | molI = currentCutoffGroup->getMember(cg); |
528 | + | |
529 | + | // tempI is atom numbering on local processor |
530 | + | tempI = molI + atomOffset; |
531 | + | |
532 | + | #ifdef IS_MPI |
533 | + | globalID = info[k].atoms[tempI]->getGlobalIndex(); |
534 | + | info[k].globalGroupMembership[globalID] = globalGroupIndex[groupOffset]; |
535 | + | #else |
536 | + | globalID = info[k].atoms[tempI]->getIndex(); |
537 | + | info[k].globalGroupMembership[globalID] = groupOffset; |
538 | + | #endif |
539 | + | myCutoffGroup->addAtom(info[k].atoms[tempI]); |
540 | + | cutoffAtomSet.insert(tempI); |
541 | + | } |
542 | + | |
543 | + | molInfo.myCutoffGroups.push_back(myCutoffGroup); |
544 | + | groupOffset++; |
545 | + | |
546 | + | }//end for (j = 0; j < molInfo.nCutoffGroups; j++) |
547 | + | |
548 | + | |
549 | + | // create a cutoff group for every atom in current molecule which |
550 | + | // does not belong to cutoffgroup defined at mdl file |
551 | + | |
552 | + | for(j = 0; j < molInfo.nAtoms; j++){ |
553 | + | |
554 | + | if(cutoffAtomSet.find(molInfo.myAtoms[j]->getIndex()) == cutoffAtomSet.end()){ |
555 | + | myCutoffGroup = new CutoffGroup(); |
556 | + | myCutoffGroup->addAtom(molInfo.myAtoms[j]); |
557 | + | |
558 | + | #ifdef IS_MPI |
559 | + | myCutoffGroup->setGlobalIndex(globalGroupIndex[groupOffset]); |
560 | + | globalID = info[k].atoms[atomOffset + j]->getGlobalIndex(); |
561 | + | info[k].globalGroupMembership[globalID] = globalGroupIndex[groupOffset]; |
562 | + | #else |
563 | + | myCutoffGroup->setGlobalIndex(groupOffset); |
564 | + | globalID = info[k].atoms[atomOffset + j]->getIndex(); |
565 | + | info[k].globalGroupMembership[globalID] = groupOffset; |
566 | + | #endif |
567 | + | molInfo.myCutoffGroups.push_back(myCutoffGroup); |
568 | + | groupOffset++; |
569 | + | } |
570 | + | } |
571 | + | |
572 | + | // After this is all set up, scan through the atoms to |
573 | + | // see if they can be added to the integrableObjects: |
574 | + | |
575 | + | molInfo.myIntegrableObjects.clear(); |
576 | + | |
577 | + | |
578 | + | for (j = 0; j < molInfo.nAtoms; j++){ |
579 | + | |
580 | + | #ifdef IS_MPI |
581 | + | slJ = molInfo.myAtoms[j]->getGlobalIndex(); |
582 | + | #else |
583 | + | slJ = j+atomOffset; |
584 | + | #endif |
585 | + | |
586 | + | // if they aren't on the skip list, then they can be integrated |
587 | + | |
588 | + | if (skipList.find(slJ) == skipList.end()) { |
589 | + | mySD = (StuntDouble *) molInfo.myAtoms[j]; |
590 | + | info[k].integrableObjects.push_back(mySD); |
591 | + | molInfo.myIntegrableObjects.push_back(mySD); |
592 | + | } |
593 | + | } |
594 | + | |
595 | + | // all rigid bodies are integrated: |
596 | + | |
597 | + | for (j = 0; j < molInfo.nRigidBodies; j++) { |
598 | + | mySD = (StuntDouble *) molInfo.myRigidBodies[j]; |
599 | + | info[k].integrableObjects.push_back(mySD); |
600 | + | molInfo.myIntegrableObjects.push_back(mySD); |
601 | + | } |
602 | + | |
603 | + | |
604 | + | /* |
605 | + | |
606 | + | //creat ConstraintPair. |
607 | + | molInfo.myConstraintPair.clear(); |
608 | + | |
609 | + | for (j = 0; j < molInfo.nBonds; j++){ |
610 | + | |
611 | + | //if both atoms are in the same rigid body, just skip it |
612 | + | currentBond = comp_stamps[stampID]->getBond(j); |
613 | + | if(!comp_stamps[stampID]->isBondInSameRigidBody(currentBond)){ |
614 | + | |
615 | + | tempI = currentBond->getA() + atomOffset; |
616 | + | if( comp_stamps[stampID]->isAtomInRigidBody(currentBond->getA(), whichRigidBody, consAtomIndex)) |
617 | + | consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex); |
618 | + | else |
619 | + | consElement1 = new ConstraintAtom(info[k].atoms[tempI]); |
620 | + | |
621 | + | tempJ = currentBond->getB() + atomOffset; |
622 | + | if(comp_stamps[stampID]->isAtomInRigidBody(currentBond->getB(), whichRigidBody, consAtomIndex)) |
623 | + | consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex); |
624 | + | else |
625 | + | consElement2 = new ConstraintAtom(info[k].atoms[tempJ]); |
626 | + | |
627 | + | consPair = new DistanceConstraintPair(consElement1, consElement2); |
628 | + | molInfo.myConstraintPairs.push_back(consPair); |
629 | + | } |
630 | + | } |
631 | + | |
632 | + | //loop over rigid bodies, if two rigid bodies share same joint, creat a HingeConstraintPair |
633 | + | for (int rb1 = 0; rb1 < molInfo.nRigidBodies -1 ; rb1++){ |
634 | + | for (int rb2 = rb1 + 1; rb2 < molInfo.nRigidBodies ; rb2++){ |
635 | + | |
636 | + | jointAtoms = comp_stamps[stampID]->getJointAtoms(rb1, rb2); |
637 | + | |
638 | + | for(size_t m = 0; m < jointAtoms.size(); m++){ |
639 | + | consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[rb1], jointAtoms[m].first); |
640 | + | consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[rb2], jointAtoms[m].second); |
641 | + | |
642 | + | consPair = new JointConstraintPair(consElement1, consElement2); |
643 | + | molInfo.myConstraintPairs.push_back(consPair); |
644 | + | } |
645 | + | |
646 | + | } |
647 | + | } |
648 | + | |
649 | + | */ |
650 | + | // send the arrays off to the forceField for init. |
651 | + | |
652 | the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms); | |
653 | the_ff->initializeBonds(molInfo.nBonds, molInfo.myBonds, theBonds); | |
654 | the_ff->initializeBends(molInfo.nBends, molInfo.myBends, theBends); | |
655 | the_ff->initializeTorsions(molInfo.nTorsions, molInfo.myTorsions, | |
656 | theTorsions); | |
657 | ||
368 | – | |
658 | info[k].molecules[i].initialize(molInfo); | |
659 | < | |
660 | < | |
659 | > | |
660 | > | |
661 | atomOffset += molInfo.nAtoms; | |
662 | delete[] theBonds; | |
663 | delete[] theBends; | |
664 | delete[] theTorsions; | |
665 | } | |
666 | + | |
667 | + | |
668 | + | |
669 | + | #ifdef IS_MPI |
670 | + | // Since the globalGroupMembership has been zero filled and we've only |
671 | + | // poked values into the atoms we know, we can do an Allreduce |
672 | + | // to get the full globalGroupMembership array (We think). |
673 | + | // This would be prettier if we could use MPI_IN_PLACE like the MPI-2 |
674 | + | // docs said we could. |
675 | + | |
676 | + | int* ggMjunk = new int[mpiSim->getNAtomsGlobal()]; |
677 | + | |
678 | + | MPI_Allreduce(info[k].globalGroupMembership, |
679 | + | ggMjunk, |
680 | + | mpiSim->getNAtomsGlobal(), |
681 | + | MPI_INT, MPI_SUM, MPI_COMM_WORLD); |
682 | + | |
683 | + | for (i = 0; i < mpiSim->getNAtomsGlobal(); i++) |
684 | + | info[k].globalGroupMembership[i] = ggMjunk[i]; |
685 | + | |
686 | + | delete[] ggMjunk; |
687 | + | |
688 | + | #endif |
689 | + | |
690 | + | |
691 | + | |
692 | } | |
693 | ||
694 | #ifdef IS_MPI | |
# | Line 381 | Line 696 | void SimSetup::makeMolecules(void){ | |
696 | MPIcheckPoint(); | |
697 | #endif // is_mpi | |
698 | ||
384 | – | // clean up the forcefield |
385 | – | |
386 | – | the_ff->calcRcut(); |
387 | – | the_ff->cleanMe(); |
699 | } | |
700 | ||
701 | void SimSetup::initFromBass(void){ | |
# | Line 585 | Line 896 | void SimSetup::gatherInfo(void){ | |
896 | else if (!strcasecmp(force_field, "EAM")){ | |
897 | ffCase = FF_EAM; | |
898 | } | |
899 | + | else if (!strcasecmp(force_field, "WATER")){ |
900 | + | ffCase = FF_H2O; |
901 | + | } |
902 | else{ | |
903 | sprintf(painCave.errMsg, "SimSetup Error. Unrecognized force field -> %s\n", | |
904 | force_field); | |
# | Line 613 | Line 927 | void SimSetup::gatherInfo(void){ | |
927 | } | |
928 | else{ | |
929 | sprintf(painCave.errMsg, | |
930 | < | "SimSetup Warning. Unrecognized Ensemble -> %s, " |
931 | < | "reverting to NVE for this simulation.\n", |
930 | > | "SimSetup Warning. Unrecognized Ensemble -> %s \n" |
931 | > | "\treverting to NVE for this simulation.\n", |
932 | ensemble); | |
933 | painCave.isFatal = 0; | |
934 | simError(); | |
# | Line 646 | Line 960 | void SimSetup::gatherInfo(void){ | |
960 | if (!the_components[i]->haveNMol()){ | |
961 | // we have a problem | |
962 | sprintf(painCave.errMsg, | |
963 | < | "SimSetup Error. No global NMol or component NMol" |
964 | < | " given. Cannot calculate the number of atoms.\n"); |
963 | > | "SimSetup Error. No global NMol or component NMol given.\n" |
964 | > | "\tCannot calculate the number of atoms.\n"); |
965 | painCave.isFatal = 1; | |
966 | simError(); | |
967 | } | |
# | Line 667 | Line 981 | void SimSetup::gatherInfo(void){ | |
981 | simError(); | |
982 | } | |
983 | ||
984 | + | //check whether sample time, status time, thermal time and reset time are divisble by dt |
985 | + | if (globals->haveSampleTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){ |
986 | + | sprintf(painCave.errMsg, |
987 | + | "Sample time is not divisible by dt.\n" |
988 | + | "\tThis will result in samples that are not uniformly\n" |
989 | + | "\tdistributed in time. If this is a problem, change\n" |
990 | + | "\tyour sampleTime variable.\n"); |
991 | + | painCave.isFatal = 0; |
992 | + | simError(); |
993 | + | } |
994 | + | |
995 | + | if (globals->haveStatusTime() && !isDivisible(globals->getStatusTime(), globals->getDt())){ |
996 | + | sprintf(painCave.errMsg, |
997 | + | "Status time is not divisible by dt.\n" |
998 | + | "\tThis will result in status reports that are not uniformly\n" |
999 | + | "\tdistributed in time. If this is a problem, change \n" |
1000 | + | "\tyour statusTime variable.\n"); |
1001 | + | painCave.isFatal = 0; |
1002 | + | simError(); |
1003 | + | } |
1004 | + | |
1005 | + | if (globals->haveThermalTime() && !isDivisible(globals->getThermalTime(), globals->getDt())){ |
1006 | + | sprintf(painCave.errMsg, |
1007 | + | "Thermal time is not divisible by dt.\n" |
1008 | + | "\tThis will result in thermalizations that are not uniformly\n" |
1009 | + | "\tdistributed in time. If this is a problem, change \n" |
1010 | + | "\tyour thermalTime variable.\n"); |
1011 | + | painCave.isFatal = 0; |
1012 | + | simError(); |
1013 | + | } |
1014 | + | |
1015 | + | if (globals->haveResetTime() && !isDivisible(globals->getResetTime(), globals->getDt())){ |
1016 | + | sprintf(painCave.errMsg, |
1017 | + | "Reset time is not divisible by dt.\n" |
1018 | + | "\tThis will result in integrator resets that are not uniformly\n" |
1019 | + | "\tdistributed in time. If this is a problem, change\n" |
1020 | + | "\tyour resetTime variable.\n"); |
1021 | + | painCave.isFatal = 0; |
1022 | + | simError(); |
1023 | + | } |
1024 | + | |
1025 | // set the status, sample, and thermal kick times | |
1026 | ||
1027 | for (i = 0; i < nInfo; i++){ | |
1028 | if (globals->haveSampleTime()){ | |
1029 | info[i].sampleTime = globals->getSampleTime(); | |
1030 | info[i].statusTime = info[i].sampleTime; | |
676 | – | info[i].thermalTime = info[i].sampleTime; |
1031 | } | |
1032 | else{ | |
1033 | info[i].sampleTime = globals->getRunTime(); | |
1034 | info[i].statusTime = info[i].sampleTime; | |
681 | – | info[i].thermalTime = info[i].sampleTime; |
1035 | } | |
1036 | ||
1037 | if (globals->haveStatusTime()){ | |
# | Line 687 | Line 1040 | void SimSetup::gatherInfo(void){ | |
1040 | ||
1041 | if (globals->haveThermalTime()){ | |
1042 | info[i].thermalTime = globals->getThermalTime(); | |
1043 | + | } else { |
1044 | + | info[i].thermalTime = globals->getRunTime(); |
1045 | } | |
1046 | ||
1047 | info[i].resetIntegrator = 0; | |
# | Line 704 | Line 1059 | void SimSetup::gatherInfo(void){ | |
1059 | ||
1060 | info[i].useInitXSstate = globals->getUseInitXSstate(); | |
1061 | info[i].orthoTolerance = globals->getOrthoBoxTolerance(); | |
1062 | < | |
1062 | > | |
1063 | > | // check for thermodynamic integration |
1064 | > | if (globals->getUseSolidThermInt() && !globals->getUseLiquidThermInt()) { |
1065 | > | if (globals->haveThermIntLambda() && globals->haveThermIntK()) { |
1066 | > | info[i].useSolidThermInt = globals->getUseSolidThermInt(); |
1067 | > | info[i].thermIntLambda = globals->getThermIntLambda(); |
1068 | > | info[i].thermIntK = globals->getThermIntK(); |
1069 | > | |
1070 | > | Restraints *myRestraint = new Restraints(tot_nmol, info[i].thermIntLambda, info[i].thermIntK); |
1071 | > | info[i].restraint = myRestraint; |
1072 | > | } |
1073 | > | else { |
1074 | > | sprintf(painCave.errMsg, |
1075 | > | "SimSetup Error:\n" |
1076 | > | "\tKeyword useSolidThermInt was set to 'true' but\n" |
1077 | > | "\tthermodynamicIntegrationLambda (and/or\n" |
1078 | > | "\tthermodynamicIntegrationK) was not specified.\n" |
1079 | > | "\tPlease provide a lambda value and k value in your .bass file.\n"); |
1080 | > | painCave.isFatal = 1; |
1081 | > | simError(); |
1082 | > | } |
1083 | > | } |
1084 | > | else if(globals->getUseLiquidThermInt()) { |
1085 | > | if (globals->getUseSolidThermInt()) { |
1086 | > | sprintf( painCave.errMsg, |
1087 | > | "SimSetup Warning: It appears that you have both solid and\n" |
1088 | > | "\tliquid thermodynamic integration activated in your .bass\n" |
1089 | > | "\tfile. To avoid confusion, specify only one technique in\n" |
1090 | > | "\tyour .bass file. Liquid-state thermodynamic integration\n" |
1091 | > | "\twill be assumed for the current simulation. If this is not\n" |
1092 | > | "\twhat you desire, set useSolidThermInt to 'true' and\n" |
1093 | > | "\tuseLiquidThermInt to 'false' in your .bass file.\n"); |
1094 | > | painCave.isFatal = 0; |
1095 | > | simError(); |
1096 | > | } |
1097 | > | if (globals->haveThermIntLambda() && globals->haveThermIntK()) { |
1098 | > | info[i].useLiquidThermInt = globals->getUseLiquidThermInt(); |
1099 | > | info[i].thermIntLambda = globals->getThermIntLambda(); |
1100 | > | info[i].thermIntK = globals->getThermIntK(); |
1101 | > | } |
1102 | > | else { |
1103 | > | sprintf(painCave.errMsg, |
1104 | > | "SimSetup Error:\n" |
1105 | > | "\tKeyword useLiquidThermInt was set to 'true' but\n" |
1106 | > | "\tthermodynamicIntegrationLambda (and/or\n" |
1107 | > | "\tthermodynamicIntegrationK) was not specified.\n" |
1108 | > | "\tPlease provide a lambda value and k value in your .bass file.\n"); |
1109 | > | painCave.isFatal = 1; |
1110 | > | simError(); |
1111 | > | } |
1112 | > | } |
1113 | > | else if(globals->haveThermIntLambda() || globals->haveThermIntK()){ |
1114 | > | sprintf(painCave.errMsg, |
1115 | > | "SimSetup Warning: If you want to use Thermodynamic\n" |
1116 | > | "\tIntegration, set useSolidThermInt or useLiquidThermInt to\n" |
1117 | > | "\t'true' in your .bass file. These keywords are set to\n" |
1118 | > | "\t'false' by default, so your lambda and/or k values are\n" |
1119 | > | "\tbeing ignored.\n"); |
1120 | > | painCave.isFatal = 0; |
1121 | > | simError(); |
1122 | > | } |
1123 | } | |
1124 | ||
1125 | //setup seed for random number generator | |
# | Line 746 | Line 1161 | void SimSetup::gatherInfo(void){ | |
1161 | for (int i = 0; i < nInfo; i++){ | |
1162 | info[i].setSeed(seedValue); | |
1163 | } | |
1164 | < | |
1164 | > | |
1165 | #ifdef IS_MPI | |
1166 | < | strcpy(checkPointMsg, "Succesfully gathered all information from Bass\n"); |
1166 | > | strcpy(checkPointMsg, "Successfully gathered all information from Bass\n"); |
1167 | MPIcheckPoint(); | |
1168 | #endif // is_mpi | |
1169 | } | |
# | Line 757 | Line 1172 | void SimSetup::finalInfoCheck(void){ | |
1172 | void SimSetup::finalInfoCheck(void){ | |
1173 | int index; | |
1174 | int usesDipoles; | |
1175 | + | int usesCharges; |
1176 | int i; | |
1177 | ||
1178 | for (i = 0; i < nInfo; i++){ | |
# | Line 768 | Line 1184 | void SimSetup::finalInfoCheck(void){ | |
1184 | usesDipoles = (info[i].atoms[index])->hasDipole(); | |
1185 | index++; | |
1186 | } | |
1187 | < | |
1187 | > | index = 0; |
1188 | > | usesCharges = 0; |
1189 | > | while ((index < info[i].n_atoms) && !usesCharges){ |
1190 | > | usesCharges= (info[i].atoms[index])->hasCharge(); |
1191 | > | index++; |
1192 | > | } |
1193 | #ifdef IS_MPI | |
1194 | int myUse = usesDipoles; | |
1195 | MPI_Allreduce(&myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
1196 | #endif //is_mpi | |
1197 | ||
1198 | < | double theEcr, theEst; |
1198 | > | double theRcut, theRsw; |
1199 | ||
1200 | + | if (globals->haveRcut()) { |
1201 | + | theRcut = globals->getRcut(); |
1202 | + | |
1203 | + | if (globals->haveRsw()) |
1204 | + | theRsw = globals->getRsw(); |
1205 | + | else |
1206 | + | theRsw = theRcut; |
1207 | + | |
1208 | + | info[i].setDefaultRcut(theRcut, theRsw); |
1209 | + | |
1210 | + | } else { |
1211 | + | |
1212 | + | the_ff->calcRcut(); |
1213 | + | theRcut = info[i].getRcut(); |
1214 | + | |
1215 | + | if (globals->haveRsw()) |
1216 | + | theRsw = globals->getRsw(); |
1217 | + | else |
1218 | + | theRsw = theRcut; |
1219 | + | |
1220 | + | info[i].setDefaultRcut(theRcut, theRsw); |
1221 | + | } |
1222 | + | |
1223 | if (globals->getUseRF()){ | |
1224 | info[i].useReactionField = 1; | |
1225 | < | |
1226 | < | if (!globals->haveECR()){ |
1225 | > | |
1226 | > | if (!globals->haveRcut()){ |
1227 | sprintf(painCave.errMsg, | |
1228 | < | "SimSetup Warning: using default value of 1/2 the smallest " |
1229 | < | "box length for the electrostaticCutoffRadius.\n" |
1230 | < | "I hope you have a very fast processor!\n"); |
1228 | > | "SimSetup Warning: No value was set for the cutoffRadius.\n" |
1229 | > | "\tOOPSE will use a default value of 15.0 angstroms" |
1230 | > | "\tfor the cutoffRadius.\n"); |
1231 | painCave.isFatal = 0; | |
1232 | simError(); | |
1233 | < | double smallest; |
790 | < | smallest = info[i].boxL[0]; |
791 | < | if (info[i].boxL[1] <= smallest) |
792 | < | smallest = info[i].boxL[1]; |
793 | < | if (info[i].boxL[2] <= smallest) |
794 | < | smallest = info[i].boxL[2]; |
795 | < | theEcr = 0.5 * smallest; |
1233 | > | theRcut = 15.0; |
1234 | } | |
1235 | else{ | |
1236 | < | theEcr = globals->getECR(); |
1236 | > | theRcut = globals->getRcut(); |
1237 | } | |
1238 | ||
1239 | < | if (!globals->haveEST()){ |
1239 | > | if (!globals->haveRsw()){ |
1240 | sprintf(painCave.errMsg, | |
1241 | < | "SimSetup Warning: using default value of 0.05 * the " |
1242 | < | "electrostaticCutoffRadius for the electrostaticSkinThickness\n"); |
1241 | > | "SimSetup Warning: No value was set for switchingRadius.\n" |
1242 | > | "\tOOPSE will use a default value of\n" |
1243 | > | "\t0.95 * cutoffRadius for the switchingRadius\n"); |
1244 | painCave.isFatal = 0; | |
1245 | simError(); | |
1246 | < | theEst = 0.05 * theEcr; |
1246 | > | theRsw = 0.95 * theRcut; |
1247 | } | |
1248 | else{ | |
1249 | < | theEst = globals->getEST(); |
1249 | > | theRsw = globals->getRsw(); |
1250 | } | |
1251 | ||
1252 | < | info[i].setDefaultEcr(theEcr, theEst); |
1252 | > | info[i].setDefaultRcut(theRcut, theRsw); |
1253 | ||
1254 | if (!globals->haveDielectric()){ | |
1255 | sprintf(painCave.errMsg, | |
1256 | < | "SimSetup Error: You are trying to use Reaction Field without" |
1257 | < | "setting a dielectric constant!\n"); |
1256 | > | "SimSetup Error: No Dielectric constant was set.\n" |
1257 | > | "\tYou are trying to use Reaction Field without" |
1258 | > | "\tsetting a dielectric constant!\n"); |
1259 | painCave.isFatal = 1; | |
1260 | simError(); | |
1261 | } | |
1262 | info[i].dielectric = globals->getDielectric(); | |
1263 | } | |
1264 | else{ | |
1265 | < | if (usesDipoles){ |
1266 | < | if (!globals->haveECR()){ |
1265 | > | if (usesDipoles || usesCharges){ |
1266 | > | |
1267 | > | if (!globals->haveRcut()){ |
1268 | sprintf(painCave.errMsg, | |
1269 | < | "SimSetup Warning: using default value of 1/2 the smallest " |
1270 | < | "box length for the electrostaticCutoffRadius.\n" |
1271 | < | "I hope you have a very fast processor!\n"); |
1269 | > | "SimSetup Warning: No value was set for the cutoffRadius.\n" |
1270 | > | "\tOOPSE will use a default value of 15.0 angstroms" |
1271 | > | "\tfor the cutoffRadius.\n"); |
1272 | painCave.isFatal = 0; | |
1273 | simError(); | |
1274 | < | double smallest; |
1275 | < | smallest = info[i].boxL[0]; |
835 | < | if (info[i].boxL[1] <= smallest) |
836 | < | smallest = info[i].boxL[1]; |
837 | < | if (info[i].boxL[2] <= smallest) |
838 | < | smallest = info[i].boxL[2]; |
839 | < | theEcr = 0.5 * smallest; |
840 | < | } |
1274 | > | theRcut = 15.0; |
1275 | > | } |
1276 | else{ | |
1277 | < | theEcr = globals->getECR(); |
1277 | > | theRcut = globals->getRcut(); |
1278 | } | |
1279 | < | |
1280 | < | if (!globals->haveEST()){ |
1279 | > | |
1280 | > | if (!globals->haveRsw()){ |
1281 | sprintf(painCave.errMsg, | |
1282 | < | "SimSetup Warning: using default value of 0.05 * the " |
1283 | < | "electrostaticCutoffRadius for the " |
1284 | < | "electrostaticSkinThickness\n"); |
1282 | > | "SimSetup Warning: No value was set for switchingRadius.\n" |
1283 | > | "\tOOPSE will use a default value of\n" |
1284 | > | "\t0.95 * cutoffRadius for the switchingRadius\n"); |
1285 | painCave.isFatal = 0; | |
1286 | simError(); | |
1287 | < | theEst = 0.05 * theEcr; |
1287 | > | theRsw = 0.95 * theRcut; |
1288 | } | |
1289 | else{ | |
1290 | < | theEst = globals->getEST(); |
1290 | > | theRsw = globals->getRsw(); |
1291 | } | |
1292 | < | |
1293 | < | info[i].setDefaultEcr(theEcr, theEst); |
1292 | > | |
1293 | > | info[i].setDefaultRcut(theRcut, theRsw); |
1294 | > | |
1295 | } | |
1296 | } | |
861 | – | |
862 | – | if( !initSuspend ) |
863 | – | info[i].checkCutOffs(); |
1297 | } | |
865 | – | |
1298 | #ifdef IS_MPI | |
1299 | strcpy(checkPointMsg, "post processing checks out"); | |
1300 | MPIcheckPoint(); | |
1301 | #endif // is_mpi | |
870 | – | } |
1302 | ||
1303 | + | // clean up the forcefield |
1304 | + | the_ff->cleanMe(); |
1305 | + | } |
1306 | + | |
1307 | void SimSetup::initSystemCoords(void){ | |
1308 | int i; | |
1309 | ||
# | Line 896 | Line 1331 | void SimSetup::initSystemCoords(void){ | |
1331 | delete fileInit; | |
1332 | } | |
1333 | else{ | |
1334 | < | #ifdef IS_MPI |
900 | < | |
1334 | > | |
1335 | // no init from bass | |
1336 | < | |
1336 | > | |
1337 | sprintf(painCave.errMsg, | |
1338 | < | "Cannot intialize a parallel simulation without an initial configuration file.\n"); |
1338 | > | "Cannot intialize a simulation without an initial configuration file.\n"); |
1339 | painCave.isFatal = 1;; | |
1340 | simError(); | |
1341 | < | |
908 | < | #else |
909 | < | |
910 | < | initFromBass(); |
911 | < | |
912 | < | |
913 | < | #endif |
1341 | > | |
1342 | } | |
1343 | ||
1344 | #ifdef IS_MPI | |
# | Line 1004 | Line 1432 | void SimSetup::makeOutNames(void){ | |
1432 | } | |
1433 | } | |
1434 | ||
1435 | + | strcpy(info[k].rawPotName, inFileName); |
1436 | + | nameLength = strlen(info[k].rawPotName); |
1437 | + | endTest = &(info[k].rawPotName[nameLength - 5]); |
1438 | + | if (!strcmp(endTest, ".bass")){ |
1439 | + | strcpy(endTest, ".raw"); |
1440 | + | } |
1441 | + | else if (!strcmp(endTest, ".BASS")){ |
1442 | + | strcpy(endTest, ".raw"); |
1443 | + | } |
1444 | + | else{ |
1445 | + | endTest = &(info[k].rawPotName[nameLength - 4]); |
1446 | + | if (!strcmp(endTest, ".bss")){ |
1447 | + | strcpy(endTest, ".raw"); |
1448 | + | } |
1449 | + | else if (!strcmp(endTest, ".mdl")){ |
1450 | + | strcpy(endTest, ".raw"); |
1451 | + | } |
1452 | + | else{ |
1453 | + | strcat(info[k].rawPotName, ".raw"); |
1454 | + | } |
1455 | + | } |
1456 | + | |
1457 | #ifdef IS_MPI | |
1458 | ||
1459 | } | |
# | Line 1064 | Line 1514 | void SimSetup::createFF(void){ | |
1514 | the_ff = new EAM_FF(); | |
1515 | break; | |
1516 | ||
1517 | + | case FF_H2O: |
1518 | + | the_ff = new WATER(); |
1519 | + | break; |
1520 | + | |
1521 | default: | |
1522 | sprintf(painCave.errMsg, | |
1523 | "SimSetup Error. Unrecognized force field in case statement.\n"); | |
# | Line 1084 | Line 1538 | void SimSetup::compList(void){ | |
1538 | LinkedMolStamp* headStamp = new LinkedMolStamp(); | |
1539 | LinkedMolStamp* currentStamp = NULL; | |
1540 | comp_stamps = new MoleculeStamp * [n_components]; | |
1541 | + | bool haveCutoffGroups; |
1542 | ||
1543 | + | haveCutoffGroups = false; |
1544 | + | |
1545 | // make an array of molecule stamps that match the components used. | |
1546 | // also extract the used stamps out into a separate linked list | |
1547 | ||
# | Line 1119 | Line 1576 | void SimSetup::compList(void){ | |
1576 | headStamp->add(currentStamp); | |
1577 | comp_stamps[i] = headStamp->match(id); | |
1578 | } | |
1579 | + | |
1580 | + | if(comp_stamps[i]->getNCutoffGroups() > 0) |
1581 | + | haveCutoffGroups = true; |
1582 | } | |
1583 | + | |
1584 | + | for (i = 0; i < nInfo; i++) |
1585 | + | info[i].haveCutoffGroups = haveCutoffGroups; |
1586 | ||
1587 | #ifdef IS_MPI | |
1588 | strcpy(checkPointMsg, "Component stamps successfully extracted\n"); | |
# | Line 1128 | Line 1591 | void SimSetup::calcSysValues(void){ | |
1591 | } | |
1592 | ||
1593 | void SimSetup::calcSysValues(void){ | |
1594 | < | int i; |
1594 | > | int i, j; |
1595 | > | int ncutgroups, atomsingroups, ngroupsinstamp; |
1596 | ||
1597 | int* molMembershipArray; | |
1598 | + | CutoffGroupStamp* cg; |
1599 | ||
1600 | tot_atoms = 0; | |
1601 | tot_bonds = 0; | |
1602 | tot_bends = 0; | |
1603 | tot_torsions = 0; | |
1604 | + | tot_rigid = 0; |
1605 | + | tot_groups = 0; |
1606 | for (i = 0; i < n_components; i++){ | |
1607 | tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms(); | |
1608 | tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds(); | |
1609 | tot_bends += components_nmol[i] * comp_stamps[i]->getNBends(); | |
1610 | tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions(); | |
1611 | < | } |
1611 | > | tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies(); |
1612 | ||
1613 | + | ncutgroups = comp_stamps[i]->getNCutoffGroups(); |
1614 | + | atomsingroups = 0; |
1615 | + | for (j=0; j < ncutgroups; j++) { |
1616 | + | cg = comp_stamps[i]->getCutoffGroup(j); |
1617 | + | atomsingroups += cg->getNMembers(); |
1618 | + | } |
1619 | + | ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups + ncutgroups; |
1620 | + | tot_groups += components_nmol[i] * ngroupsinstamp; |
1621 | + | } |
1622 | + | |
1623 | tot_SRI = tot_bonds + tot_bends + tot_torsions; | |
1624 | molMembershipArray = new int[tot_atoms]; | |
1625 | ||
# | Line 1153 | Line 1630 | void SimSetup::calcSysValues(void){ | |
1630 | info[i].n_torsions = tot_torsions; | |
1631 | info[i].n_SRI = tot_SRI; | |
1632 | info[i].n_mol = tot_nmol; | |
1633 | < | |
1633 | > | info[i].ngroup = tot_groups; |
1634 | info[i].molMembershipArray = molMembershipArray; | |
1635 | } | |
1636 | } | |
# | Line 1164 | Line 1641 | void SimSetup::mpiMolDivide(void){ | |
1641 | int i, j, k; | |
1642 | int localMol, allMol; | |
1643 | int local_atoms, local_bonds, local_bends, local_torsions, local_SRI; | |
1644 | + | int local_rigid, local_groups; |
1645 | + | vector<int> globalMolIndex; |
1646 | + | int ncutgroups, atomsingroups, ngroupsinstamp; |
1647 | + | CutoffGroupStamp* cg; |
1648 | ||
1649 | mpiSim = new mpiSimulation(info); | |
1650 | ||
1651 | < | globalIndex = mpiSim->divideLabor(); |
1651 | > | mpiSim->divideLabor(); |
1652 | > | globalAtomIndex = mpiSim->getGlobalAtomIndex(); |
1653 | > | globalGroupIndex = mpiSim->getGlobalGroupIndex(); |
1654 | > | //globalMolIndex = mpiSim->getGlobalMolIndex(); |
1655 | ||
1656 | // set up the local variables | |
1657 | ||
# | Line 1180 | Line 1664 | void SimSetup::mpiMolDivide(void){ | |
1664 | local_bonds = 0; | |
1665 | local_bends = 0; | |
1666 | local_torsions = 0; | |
1667 | < | globalAtomIndex = 0; |
1667 | > | local_rigid = 0; |
1668 | > | local_groups = 0; |
1669 | > | globalAtomCounter = 0; |
1670 | ||
1185 | – | |
1671 | for (i = 0; i < n_components; i++){ | |
1672 | for (j = 0; j < components_nmol[i]; j++){ | |
1673 | if (mol2proc[allMol] == worldRank){ | |
# | Line 1190 | Line 1675 | void SimSetup::mpiMolDivide(void){ | |
1675 | local_bonds += comp_stamps[i]->getNBonds(); | |
1676 | local_bends += comp_stamps[i]->getNBends(); | |
1677 | local_torsions += comp_stamps[i]->getNTorsions(); | |
1678 | + | local_rigid += comp_stamps[i]->getNRigidBodies(); |
1679 | + | |
1680 | + | ncutgroups = comp_stamps[i]->getNCutoffGroups(); |
1681 | + | atomsingroups = 0; |
1682 | + | for (k=0; k < ncutgroups; k++) { |
1683 | + | cg = comp_stamps[i]->getCutoffGroup(k); |
1684 | + | atomsingroups += cg->getNMembers(); |
1685 | + | } |
1686 | + | ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups + |
1687 | + | ncutgroups; |
1688 | + | local_groups += ngroupsinstamp; |
1689 | + | |
1690 | localMol++; | |
1691 | } | |
1692 | for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){ | |
1693 | < | info[0].molMembershipArray[globalAtomIndex] = allMol; |
1694 | < | globalAtomIndex++; |
1693 | > | info[0].molMembershipArray[globalAtomCounter] = allMol; |
1694 | > | globalAtomCounter++; |
1695 | } | |
1696 | ||
1697 | allMol++; | |
# | Line 1202 | Line 1699 | void SimSetup::mpiMolDivide(void){ | |
1699 | } | |
1700 | local_SRI = local_bonds + local_bends + local_torsions; | |
1701 | ||
1702 | < | info[0].n_atoms = mpiSim->getMyNlocal(); |
1703 | < | |
1702 | > | info[0].n_atoms = mpiSim->getNAtomsLocal(); |
1703 | > | |
1704 | if (local_atoms != info[0].n_atoms){ | |
1705 | sprintf(painCave.errMsg, | |
1706 | < | "SimSetup error: mpiSim's localAtom (%d) and SimSetup's" |
1707 | < | " localAtom (%d) are not equal.\n", |
1706 | > | "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n" |
1707 | > | "\tlocalAtom (%d) are not equal.\n", |
1708 | info[0].n_atoms, local_atoms); | |
1709 | painCave.isFatal = 1; | |
1710 | simError(); | |
1711 | } | |
1712 | ||
1713 | + | info[0].ngroup = mpiSim->getNGroupsLocal(); |
1714 | + | if (local_groups != info[0].ngroup){ |
1715 | + | sprintf(painCave.errMsg, |
1716 | + | "SimSetup error: mpiSim's localGroups (%d) and SimSetup's\n" |
1717 | + | "\tlocalGroups (%d) are not equal.\n", |
1718 | + | info[0].ngroup, local_groups); |
1719 | + | painCave.isFatal = 1; |
1720 | + | simError(); |
1721 | + | } |
1722 | + | |
1723 | info[0].n_bonds = local_bonds; | |
1724 | info[0].n_bends = local_bends; | |
1725 | info[0].n_torsions = local_torsions; | |
# | Line 1235 | Line 1742 | void SimSetup::makeSysArrays(void){ | |
1742 | ||
1743 | Atom** the_atoms; | |
1744 | Molecule* the_molecules; | |
1238 | – | Exclude** the_excludes; |
1745 | ||
1240 | – | |
1746 | for (l = 0; l < nInfo; l++){ | |
1747 | // create the atom and short range interaction arrays | |
1748 | ||
# | Line 1251 | Line 1756 | void SimSetup::makeSysArrays(void){ | |
1756 | ||
1757 | ||
1758 | molIndex = 0; | |
1759 | < | for (i = 0; i < mpiSim->getTotNmol(); i++){ |
1759 | > | for (i = 0; i < mpiSim->getNMolGlobal(); i++){ |
1760 | if (mol2proc[i] == worldRank){ | |
1761 | the_molecules[molIndex].setStampID(molCompType[i]); | |
1762 | the_molecules[molIndex].setMyIndex(molIndex); | |
# | Line 1263 | Line 1768 | void SimSetup::makeSysArrays(void){ | |
1768 | #else // is_mpi | |
1769 | ||
1770 | molIndex = 0; | |
1771 | < | globalAtomIndex = 0; |
1771 | > | globalAtomCounter = 0; |
1772 | for (i = 0; i < n_components; i++){ | |
1773 | for (j = 0; j < components_nmol[i]; j++){ | |
1774 | the_molecules[molIndex].setStampID(i); | |
1775 | the_molecules[molIndex].setMyIndex(molIndex); | |
1776 | the_molecules[molIndex].setGlobalIndex(molIndex); | |
1777 | for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){ | |
1778 | < | info[l].molMembershipArray[globalAtomIndex] = molIndex; |
1779 | < | globalAtomIndex++; |
1778 | > | info[l].molMembershipArray[globalAtomCounter] = molIndex; |
1779 | > | globalAtomCounter++; |
1780 | } | |
1781 | molIndex++; | |
1782 | } | |
# | Line 1280 | Line 1785 | void SimSetup::makeSysArrays(void){ | |
1785 | ||
1786 | #endif // is_mpi | |
1787 | ||
1788 | < | |
1789 | < | if (info[l].n_SRI){ |
1790 | < | Exclude::createArray(info[l].n_SRI); |
1286 | < | the_excludes = new Exclude * [info[l].n_SRI]; |
1287 | < | for (int ex = 0; ex < info[l].n_SRI; ex++){ |
1288 | < | the_excludes[ex] = new Exclude(ex); |
1289 | < | } |
1290 | < | info[l].globalExcludes = new int; |
1291 | < | info[l].n_exclude = info[l].n_SRI; |
1292 | < | } |
1293 | < | else{ |
1294 | < | Exclude::createArray(1); |
1295 | < | the_excludes = new Exclude * ; |
1296 | < | the_excludes[0] = new Exclude(0); |
1297 | < | the_excludes[0]->setPair(0, 0); |
1298 | < | info[l].globalExcludes = new int; |
1299 | < | info[l].globalExcludes[0] = 0; |
1300 | < | info[l].n_exclude = 0; |
1301 | < | } |
1302 | < | |
1788 | > | info[l].globalExcludes = new int; |
1789 | > | info[l].globalExcludes[0] = 0; |
1790 | > | |
1791 | // set the arrays into the SimInfo object | |
1792 | ||
1793 | info[l].atoms = the_atoms; | |
1794 | info[l].molecules = the_molecules; | |
1795 | info[l].nGlobalExcludes = 0; | |
1796 | < | info[l].excludes = the_excludes; |
1309 | < | |
1796 | > | |
1797 | the_ff->setSimInfo(info); | |
1798 | } | |
1799 | } | |
# | Line 1349 | Line 1836 | void SimSetup::makeIntegrator(void){ | |
1836 | else{ | |
1837 | sprintf(painCave.errMsg, | |
1838 | "SimSetup error: If you use the NVT\n" | |
1839 | < | " ensemble, you must set tauThermostat.\n"); |
1839 | > | "\tensemble, you must set tauThermostat.\n"); |
1840 | painCave.isFatal = 1; | |
1841 | simError(); | |
1842 | } | |
# | Line 1372 | Line 1859 | void SimSetup::makeIntegrator(void){ | |
1859 | else{ | |
1860 | sprintf(painCave.errMsg, | |
1861 | "SimSetup error: If you use a constant pressure\n" | |
1862 | < | " ensemble, you must set targetPressure in the BASS file.\n"); |
1862 | > | "\tensemble, you must set targetPressure in the BASS file.\n"); |
1863 | painCave.isFatal = 1; | |
1864 | simError(); | |
1865 | } | |
# | Line 1382 | Line 1869 | void SimSetup::makeIntegrator(void){ | |
1869 | else{ | |
1870 | sprintf(painCave.errMsg, | |
1871 | "SimSetup error: If you use an NPT\n" | |
1872 | < | " ensemble, you must set tauThermostat.\n"); |
1872 | > | "\tensemble, you must set tauThermostat.\n"); |
1873 | painCave.isFatal = 1; | |
1874 | simError(); | |
1875 | } | |
# | Line 1392 | Line 1879 | void SimSetup::makeIntegrator(void){ | |
1879 | else{ | |
1880 | sprintf(painCave.errMsg, | |
1881 | "SimSetup error: If you use an NPT\n" | |
1882 | < | " ensemble, you must set tauBarostat.\n"); |
1882 | > | "\tensemble, you must set tauBarostat.\n"); |
1883 | painCave.isFatal = 1; | |
1884 | simError(); | |
1885 | } | |
# | Line 1415 | Line 1902 | void SimSetup::makeIntegrator(void){ | |
1902 | else{ | |
1903 | sprintf(painCave.errMsg, | |
1904 | "SimSetup error: If you use a constant pressure\n" | |
1905 | < | " ensemble, you must set targetPressure in the BASS file.\n"); |
1905 | > | "\tensemble, you must set targetPressure in the BASS file.\n"); |
1906 | painCave.isFatal = 1; | |
1907 | simError(); | |
1908 | } | |
# | Line 1426 | Line 1913 | void SimSetup::makeIntegrator(void){ | |
1913 | else{ | |
1914 | sprintf(painCave.errMsg, | |
1915 | "SimSetup error: If you use an NPT\n" | |
1916 | < | " ensemble, you must set tauThermostat.\n"); |
1916 | > | "\tensemble, you must set tauThermostat.\n"); |
1917 | painCave.isFatal = 1; | |
1918 | simError(); | |
1919 | } | |
# | Line 1437 | Line 1924 | void SimSetup::makeIntegrator(void){ | |
1924 | else{ | |
1925 | sprintf(painCave.errMsg, | |
1926 | "SimSetup error: If you use an NPT\n" | |
1927 | < | " ensemble, you must set tauBarostat.\n"); |
1927 | > | "\tensemble, you must set tauBarostat.\n"); |
1928 | painCave.isFatal = 1; | |
1929 | simError(); | |
1930 | } | |
# | Line 1460 | Line 1947 | void SimSetup::makeIntegrator(void){ | |
1947 | else{ | |
1948 | sprintf(painCave.errMsg, | |
1949 | "SimSetup error: If you use a constant pressure\n" | |
1950 | < | " ensemble, you must set targetPressure in the BASS file.\n"); |
1950 | > | "\tensemble, you must set targetPressure in the BASS file.\n"); |
1951 | painCave.isFatal = 1; | |
1952 | simError(); | |
1953 | } | |
# | Line 1470 | Line 1957 | void SimSetup::makeIntegrator(void){ | |
1957 | else{ | |
1958 | sprintf(painCave.errMsg, | |
1959 | "SimSetup error: If you use an NPT\n" | |
1960 | < | " ensemble, you must set tauThermostat.\n"); |
1960 | > | "\tensemble, you must set tauThermostat.\n"); |
1961 | painCave.isFatal = 1; | |
1962 | simError(); | |
1963 | } | |
# | Line 1480 | Line 1967 | void SimSetup::makeIntegrator(void){ | |
1967 | else{ | |
1968 | sprintf(painCave.errMsg, | |
1969 | "SimSetup error: If you use an NPT\n" | |
1970 | < | " ensemble, you must set tauBarostat.\n"); |
1970 | > | "\tensemble, you must set tauBarostat.\n"); |
1971 | painCave.isFatal = 1; | |
1972 | simError(); | |
1973 | } | |
# | Line 1533 | Line 2020 | void SimSetup::setupZConstraint(SimInfo& theInfo){ | |
2020 | } | |
2021 | else{ | |
2022 | sprintf(painCave.errMsg, | |
2023 | < | "ZConstraint error: If you use an ZConstraint\n" |
2024 | < | " , you must set sample time.\n"); |
2023 | > | "ZConstraint error: If you use a ZConstraint,\n" |
2024 | > | "\tyou must set zconsTime.\n"); |
2025 | painCave.isFatal = 1; | |
2026 | simError(); | |
2027 | } | |
# | Line 1549 | Line 2036 | void SimSetup::setupZConstraint(SimInfo& theInfo){ | |
2036 | else{ | |
2037 | double defaultZConsTol = 0.01; | |
2038 | sprintf(painCave.errMsg, | |
2039 | < | "ZConstraint Waring: Tolerance for z-constraint methodl is not specified\n" |
2040 | < | " , default value %f is used.\n", |
2039 | > | "ZConstraint Warning: Tolerance for z-constraint method is not specified.\n" |
2040 | > | "\tOOPSE will use a default value of %f.\n" |
2041 | > | "\tTo set the tolerance, use the zconsTol variable.\n", |
2042 | defaultZConsTol); | |
2043 | painCave.isFatal = 0; | |
2044 | simError(); | |
# | Line 1568 | Line 2056 | void SimSetup::setupZConstraint(SimInfo& theInfo){ | |
2056 | } | |
2057 | else{ | |
2058 | sprintf(painCave.errMsg, | |
2059 | < | "ZConstraint Warning: User does not set force Subtraction policy, " |
2060 | < | "PolicyByMass is used\n"); |
2059 | > | "ZConstraint Warning: No force subtraction policy was set.\n" |
2060 | > | "\tOOPSE will use PolicyByMass.\n" |
2061 | > | "\tTo set the policy, use the zconsForcePolicy variable.\n"); |
2062 | painCave.isFatal = 0; | |
2063 | simError(); | |
2064 | zconsForcePolicy->setData("BYMASS"); | |
# | Line 1577 | Line 2066 | void SimSetup::setupZConstraint(SimInfo& theInfo){ | |
2066 | ||
2067 | theInfo.addProperty(zconsForcePolicy); | |
2068 | ||
2069 | + | //set zcons gap |
2070 | + | DoubleData* zconsGap = new DoubleData(); |
2071 | + | zconsGap->setID(ZCONSGAP_ID); |
2072 | + | |
2073 | + | if (globals->haveZConsGap()){ |
2074 | + | zconsGap->setData(globals->getZconsGap()); |
2075 | + | theInfo.addProperty(zconsGap); |
2076 | + | } |
2077 | + | |
2078 | + | //set zcons fixtime |
2079 | + | DoubleData* zconsFixtime = new DoubleData(); |
2080 | + | zconsFixtime->setID(ZCONSFIXTIME_ID); |
2081 | + | |
2082 | + | if (globals->haveZConsFixTime()){ |
2083 | + | zconsFixtime->setData(globals->getZconsFixtime()); |
2084 | + | theInfo.addProperty(zconsFixtime); |
2085 | + | } |
2086 | + | |
2087 | + | //set zconsUsingSMD |
2088 | + | IntData* zconsUsingSMD = new IntData(); |
2089 | + | zconsUsingSMD->setID(ZCONSUSINGSMD_ID); |
2090 | + | |
2091 | + | if (globals->haveZConsUsingSMD()){ |
2092 | + | zconsUsingSMD->setData(globals->getZconsUsingSMD()); |
2093 | + | theInfo.addProperty(zconsUsingSMD); |
2094 | + | } |
2095 | + | |
2096 | //Determine the name of ouput file and add it into SimInfo's property list | |
2097 | //Be careful, do not use inFileName, since it is a pointer which | |
2098 | //point to a string at master node, and slave nodes do not contain that string | |
# | Line 1606 | Line 2122 | void SimSetup::setupZConstraint(SimInfo& theInfo){ | |
2122 | tempParaItem.zPos = zconStamp[i]->getZpos(); | |
2123 | tempParaItem.zconsIndex = zconStamp[i]->getMolIndex(); | |
2124 | tempParaItem.kRatio = zconStamp[i]->getKratio(); | |
2125 | < | |
2125 | > | tempParaItem.havingCantVel = zconStamp[i]->haveCantVel(); |
2126 | > | tempParaItem.cantVel = zconStamp[i]->getCantVel(); |
2127 | zconsParaData->addItem(tempParaItem); | |
2128 | } | |
2129 | ||
2130 | //check the uniqueness of index | |
2131 | if(!zconsParaData->isIndexUnique()){ | |
2132 | sprintf(painCave.errMsg, | |
2133 | < | "ZConstraint Error: molIndex is not unique\n"); |
2133 | > | "ZConstraint Error: molIndex is not unique!\n"); |
2134 | painCave.isFatal = 1; | |
2135 | simError(); | |
2136 | } | |
# | Line 1623 | Line 2140 | void SimSetup::setupZConstraint(SimInfo& theInfo){ | |
2140 | ||
2141 | //push data into siminfo, therefore, we can retrieve later | |
2142 | theInfo.addProperty(zconsParaData); | |
2143 | + | } |
2144 | + | |
2145 | + | void SimSetup::makeMinimizer(){ |
2146 | + | |
2147 | + | OOPSEMinimizer* myOOPSEMinimizer; |
2148 | + | MinimizerParameterSet* param; |
2149 | + | char minimizerName[100]; |
2150 | + | |
2151 | + | for (int i = 0; i < nInfo; i++){ |
2152 | + | |
2153 | + | //prepare parameter set for minimizer |
2154 | + | param = new MinimizerParameterSet(); |
2155 | + | param->setDefaultParameter(); |
2156 | + | |
2157 | + | if (globals->haveMinimizer()){ |
2158 | + | param->setFTol(globals->getMinFTol()); |
2159 | + | } |
2160 | + | |
2161 | + | if (globals->haveMinGTol()){ |
2162 | + | param->setGTol(globals->getMinGTol()); |
2163 | + | } |
2164 | + | |
2165 | + | if (globals->haveMinMaxIter()){ |
2166 | + | param->setMaxIteration(globals->getMinMaxIter()); |
2167 | + | } |
2168 | + | |
2169 | + | if (globals->haveMinWriteFrq()){ |
2170 | + | param->setMaxIteration(globals->getMinMaxIter()); |
2171 | + | } |
2172 | + | |
2173 | + | if (globals->haveMinWriteFrq()){ |
2174 | + | param->setWriteFrq(globals->getMinWriteFrq()); |
2175 | + | } |
2176 | + | |
2177 | + | if (globals->haveMinStepSize()){ |
2178 | + | param->setStepSize(globals->getMinStepSize()); |
2179 | + | } |
2180 | + | |
2181 | + | if (globals->haveMinLSMaxIter()){ |
2182 | + | param->setLineSearchMaxIteration(globals->getMinLSMaxIter()); |
2183 | + | } |
2184 | + | |
2185 | + | if (globals->haveMinLSTol()){ |
2186 | + | param->setLineSearchTol(globals->getMinLSTol()); |
2187 | + | } |
2188 | + | |
2189 | + | strcpy(minimizerName, globals->getMinimizer()); |
2190 | + | |
2191 | + | if (!strcasecmp(minimizerName, "CG")){ |
2192 | + | myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param); |
2193 | + | } |
2194 | + | else if (!strcasecmp(minimizerName, "SD")){ |
2195 | + | //myOOPSEMinimizer = MinimizerFactory.creatMinimizer("", &(info[i]), the_ff, param); |
2196 | + | myOOPSEMinimizer = new SDMinimizer(&(info[i]), the_ff, param); |
2197 | + | } |
2198 | + | else{ |
2199 | + | sprintf(painCave.errMsg, |
2200 | + | "SimSetup error: Unrecognized Minimizer, use Conjugate Gradient \n"); |
2201 | + | painCave.isFatal = 0; |
2202 | + | simError(); |
2203 | + | |
2204 | + | myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param); |
2205 | + | } |
2206 | + | info[i].the_integrator = myOOPSEMinimizer; |
2207 | + | |
2208 | + | //store the minimizer into simInfo |
2209 | + | info[i].the_minimizer = myOOPSEMinimizer; |
2210 | + | info[i].has_minimizer = true; |
2211 | + | } |
2212 | + | |
2213 | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |