9 |
|
#include "parse_me.h" |
10 |
|
#include "Integrator.hpp" |
11 |
|
#include "simError.h" |
12 |
+ |
#include "RigidBody.hpp" |
13 |
+ |
//#include "ConjugateMinimizer.hpp" |
14 |
+ |
#include "OOPSEMinimizer.hpp" |
15 |
|
|
16 |
|
#ifdef IS_MPI |
17 |
|
#include "mpiBASS.h" |
27 |
|
#define NPTxyz_ENS 4 |
28 |
|
|
29 |
|
|
30 |
< |
#define FF_DUFF 0 |
31 |
< |
#define FF_LJ 1 |
32 |
< |
#define FF_EAM 2 |
30 |
> |
#define FF_DUFF 0 |
31 |
> |
#define FF_LJ 1 |
32 |
> |
#define FF_EAM 2 |
33 |
> |
#define FF_H2O 3 |
34 |
|
|
35 |
|
using namespace std; |
36 |
|
|
148 |
|
|
149 |
|
makeOutNames(); |
150 |
|
|
151 |
< |
// make the integrator |
152 |
< |
|
153 |
< |
makeIntegrator(); |
154 |
< |
|
151 |
> |
if (globals->haveMinimizer()) |
152 |
> |
// make minimizer |
153 |
> |
makeMinimizer(); |
154 |
> |
else |
155 |
> |
// make the integrator |
156 |
> |
makeIntegrator(); |
157 |
> |
|
158 |
|
#ifdef IS_MPI |
159 |
|
mpiSim->mpiRefresh(); |
160 |
|
#endif |
166 |
|
|
167 |
|
|
168 |
|
void SimSetup::makeMolecules(void){ |
169 |
< |
int k; |
170 |
< |
int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset; |
169 |
> |
int i, j, k; |
170 |
> |
int exI, exJ, exK, exL, slI; |
171 |
> |
int tempI, tempJ, tempK, tempL; |
172 |
> |
int molI; |
173 |
> |
int stampID, atomOffset, rbOffset; |
174 |
|
molInit molInfo; |
175 |
|
DirectionalAtom* dAtom; |
176 |
+ |
RigidBody* myRB; |
177 |
+ |
StuntDouble* mySD; |
178 |
|
LinkedAssign* extras; |
179 |
|
LinkedAssign* current_extra; |
180 |
|
AtomStamp* currentAtom; |
181 |
|
BondStamp* currentBond; |
182 |
|
BendStamp* currentBend; |
183 |
|
TorsionStamp* currentTorsion; |
184 |
+ |
RigidBodyStamp* currentRigidBody; |
185 |
|
|
186 |
|
bond_pair* theBonds; |
187 |
|
bend_set* theBends; |
188 |
|
torsion_set* theTorsions; |
189 |
|
|
190 |
+ |
set<int> skipList; |
191 |
|
|
192 |
+ |
double phi, theta, psi; |
193 |
+ |
|
194 |
|
//init the forceField paramters |
195 |
|
|
196 |
|
the_ff->readParams(); |
197 |
|
|
182 |
– |
|
198 |
|
// init the atoms |
199 |
|
|
200 |
< |
double ux, uy, uz, u, uSqr; |
200 |
> |
int nMembers, nNew, rb1, rb2; |
201 |
|
|
202 |
|
for (k = 0; k < nInfo; k++){ |
203 |
|
the_ff->setSimInfo(&(info[k])); |
204 |
|
|
205 |
|
atomOffset = 0; |
206 |
< |
excludeOffset = 0; |
206 |
> |
|
207 |
|
for (i = 0; i < info[k].n_mol; i++){ |
208 |
|
stampID = info[k].molecules[i].getStampID(); |
209 |
|
|
211 |
|
molInfo.nBonds = comp_stamps[stampID]->getNBonds(); |
212 |
|
molInfo.nBends = comp_stamps[stampID]->getNBends(); |
213 |
|
molInfo.nTorsions = comp_stamps[stampID]->getNTorsions(); |
214 |
< |
molInfo.nExcludes = molInfo.nBonds + molInfo.nBends + molInfo.nTorsions; |
215 |
< |
|
214 |
> |
molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies(); |
215 |
> |
|
216 |
|
molInfo.myAtoms = &(info[k].atoms[atomOffset]); |
202 |
– |
molInfo.myExcludes = &(info[k].excludes[excludeOffset]); |
203 |
– |
molInfo.myBonds = new Bond * [molInfo.nBonds]; |
204 |
– |
molInfo.myBends = new Bend * [molInfo.nBends]; |
205 |
– |
molInfo.myTorsions = new Torsion * [molInfo.nTorsions]; |
217 |
|
|
218 |
+ |
if (molInfo.nBonds > 0) |
219 |
+ |
molInfo.myBonds = new (Bond *) [molInfo.nBonds]; |
220 |
+ |
else |
221 |
+ |
molInfo.myBonds = NULL; |
222 |
+ |
|
223 |
+ |
if (molInfo.nBends > 0) |
224 |
+ |
molInfo.myBends = new (Bend *) [molInfo.nBends]; |
225 |
+ |
else |
226 |
+ |
molInfo.myBends = NULL; |
227 |
+ |
|
228 |
+ |
if (molInfo.nTorsions > 0) |
229 |
+ |
molInfo.myTorsions = new (Torsion *) [molInfo.nTorsions]; |
230 |
+ |
else |
231 |
+ |
molInfo.myTorsions = NULL; |
232 |
+ |
|
233 |
|
theBonds = new bond_pair[molInfo.nBonds]; |
234 |
|
theBends = new bend_set[molInfo.nBends]; |
235 |
|
theTorsions = new torsion_set[molInfo.nTorsions]; |
236 |
< |
|
236 |
> |
|
237 |
|
// make the Atoms |
238 |
|
|
239 |
|
for (j = 0; j < molInfo.nAtoms; j++){ |
240 |
|
currentAtom = comp_stamps[stampID]->getAtom(j); |
241 |
+ |
|
242 |
|
if (currentAtom->haveOrientation()){ |
243 |
|
dAtom = new DirectionalAtom((j + atomOffset), |
244 |
|
info[k].getConfiguration()); |
245 |
|
info[k].n_oriented++; |
246 |
|
molInfo.myAtoms[j] = dAtom; |
247 |
|
|
248 |
< |
ux = currentAtom->getOrntX(); |
249 |
< |
uy = currentAtom->getOrntY(); |
250 |
< |
uz = currentAtom->getOrntZ(); |
248 |
> |
// Directional Atoms have standard unit vectors which are oriented |
249 |
> |
// in space using the three Euler angles. We assume the standard |
250 |
> |
// unit vector was originally along the z axis below. |
251 |
|
|
252 |
< |
uSqr = (ux * ux) + (uy * uy) + (uz * uz); |
252 |
> |
phi = currentAtom->getEulerPhi() * M_PI / 180.0; |
253 |
> |
theta = currentAtom->getEulerTheta() * M_PI / 180.0; |
254 |
> |
psi = currentAtom->getEulerPsi()* M_PI / 180.0; |
255 |
|
|
256 |
< |
u = sqrt(uSqr); |
257 |
< |
ux = ux / u; |
229 |
< |
uy = uy / u; |
230 |
< |
uz = uz / u; |
231 |
< |
|
232 |
< |
dAtom->setSUx(ux); |
233 |
< |
dAtom->setSUy(uy); |
234 |
< |
dAtom->setSUz(uz); |
256 |
> |
dAtom->setUnitFrameFromEuler(phi, theta, psi); |
257 |
> |
|
258 |
|
} |
259 |
|
else{ |
260 |
< |
molInfo.myAtoms[j] = new GeneralAtom((j + atomOffset), |
261 |
< |
info[k].getConfiguration()); |
260 |
> |
|
261 |
> |
molInfo.myAtoms[j] = new Atom((j + atomOffset), info[k].getConfiguration()); |
262 |
|
} |
263 |
+ |
|
264 |
|
molInfo.myAtoms[j]->setType(currentAtom->getType()); |
265 |
|
|
266 |
|
#ifdef IS_MPI |
276 |
|
theBonds[j].a = currentBond->getA() + atomOffset; |
277 |
|
theBonds[j].b = currentBond->getB() + atomOffset; |
278 |
|
|
279 |
< |
exI = theBonds[j].a; |
280 |
< |
exJ = theBonds[j].b; |
279 |
> |
tempI = theBonds[j].a; |
280 |
> |
tempJ = theBonds[j].b; |
281 |
|
|
258 |
– |
// exclude_I must always be the smaller of the pair |
259 |
– |
if (exI > exJ){ |
260 |
– |
tempEx = exI; |
261 |
– |
exI = exJ; |
262 |
– |
exJ = tempEx; |
263 |
– |
} |
282 |
|
#ifdef IS_MPI |
283 |
< |
tempEx = exI; |
284 |
< |
exI = info[k].atoms[tempEx]->getGlobalIndex() + 1; |
285 |
< |
tempEx = exJ; |
286 |
< |
exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1; |
283 |
> |
exI = info[k].atoms[tempI]->getGlobalIndex() + 1; |
284 |
> |
exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1; |
285 |
> |
#else |
286 |
> |
exI = tempI + 1; |
287 |
> |
exJ = tempJ + 1; |
288 |
> |
#endif |
289 |
|
|
290 |
< |
info[k].excludes[j + excludeOffset]->setPair(exI, exJ); |
271 |
< |
#else // isn't MPI |
272 |
< |
|
273 |
< |
info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1)); |
274 |
< |
#endif //is_mpi |
290 |
> |
info[k].excludes->addPair(exI, exJ); |
291 |
|
} |
276 |
– |
excludeOffset += molInfo.nBonds; |
292 |
|
|
293 |
|
//make the bends |
294 |
|
for (j = 0; j < molInfo.nBends; j++){ |
336 |
|
|
337 |
|
current_extra = current_extra->getNext(); |
338 |
|
} |
324 |
– |
} |
325 |
– |
|
326 |
– |
if (!theBends[j].isGhost){ |
327 |
– |
exI = theBends[j].a; |
328 |
– |
exJ = theBends[j].c; |
339 |
|
} |
330 |
– |
else{ |
331 |
– |
exI = theBends[j].a; |
332 |
– |
exJ = theBends[j].b; |
333 |
– |
} |
340 |
|
|
341 |
< |
// exclude_I must always be the smaller of the pair |
342 |
< |
if (exI > exJ){ |
343 |
< |
tempEx = exI; |
344 |
< |
exI = exJ; |
345 |
< |
exJ = tempEx; |
340 |
< |
} |
341 |
> |
if (theBends[j].isGhost) { |
342 |
> |
|
343 |
> |
tempI = theBends[j].a; |
344 |
> |
tempJ = theBends[j].b; |
345 |
> |
|
346 |
|
#ifdef IS_MPI |
347 |
< |
tempEx = exI; |
348 |
< |
exI = info[k].atoms[tempEx]->getGlobalIndex() + 1; |
349 |
< |
tempEx = exJ; |
350 |
< |
exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1; |
347 |
> |
exI = info[k].atoms[tempI]->getGlobalIndex() + 1; |
348 |
> |
exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1; |
349 |
> |
#else |
350 |
> |
exI = tempI + 1; |
351 |
> |
exJ = tempJ + 1; |
352 |
> |
#endif |
353 |
> |
info[k].excludes->addPair(exI, exJ); |
354 |
|
|
355 |
< |
info[k].excludes[j + excludeOffset]->setPair(exI, exJ); |
356 |
< |
#else // isn't MPI |
357 |
< |
info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1)); |
358 |
< |
#endif //is_mpi |
355 |
> |
} else { |
356 |
> |
|
357 |
> |
tempI = theBends[j].a; |
358 |
> |
tempJ = theBends[j].b; |
359 |
> |
tempK = theBends[j].c; |
360 |
> |
|
361 |
> |
#ifdef IS_MPI |
362 |
> |
exI = info[k].atoms[tempI]->getGlobalIndex() + 1; |
363 |
> |
exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1; |
364 |
> |
exK = info[k].atoms[tempK]->getGlobalIndex() + 1; |
365 |
> |
#else |
366 |
> |
exI = tempI + 1; |
367 |
> |
exJ = tempJ + 1; |
368 |
> |
exK = tempK + 1; |
369 |
> |
#endif |
370 |
> |
|
371 |
> |
info[k].excludes->addPair(exI, exK); |
372 |
> |
info[k].excludes->addPair(exI, exJ); |
373 |
> |
info[k].excludes->addPair(exJ, exK); |
374 |
> |
} |
375 |
|
} |
352 |
– |
excludeOffset += molInfo.nBends; |
376 |
|
|
377 |
|
for (j = 0; j < molInfo.nTorsions; j++){ |
378 |
|
currentTorsion = comp_stamps[stampID]->getTorsion(j); |
381 |
|
theTorsions[j].c = currentTorsion->getC() + atomOffset; |
382 |
|
theTorsions[j].d = currentTorsion->getD() + atomOffset; |
383 |
|
|
384 |
< |
exI = theTorsions[j].a; |
385 |
< |
exJ = theTorsions[j].d; |
384 |
> |
tempI = theTorsions[j].a; |
385 |
> |
tempJ = theTorsions[j].b; |
386 |
> |
tempK = theTorsions[j].c; |
387 |
> |
tempL = theTorsions[j].d; |
388 |
|
|
364 |
– |
// exclude_I must always be the smaller of the pair |
365 |
– |
if (exI > exJ){ |
366 |
– |
tempEx = exI; |
367 |
– |
exI = exJ; |
368 |
– |
exJ = tempEx; |
369 |
– |
} |
389 |
|
#ifdef IS_MPI |
390 |
< |
tempEx = exI; |
391 |
< |
exI = info[k].atoms[tempEx]->getGlobalIndex() + 1; |
392 |
< |
tempEx = exJ; |
393 |
< |
exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1; |
390 |
> |
exI = info[k].atoms[tempI]->getGlobalIndex() + 1; |
391 |
> |
exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1; |
392 |
> |
exK = info[k].atoms[tempK]->getGlobalIndex() + 1; |
393 |
> |
exL = info[k].atoms[tempL]->getGlobalIndex() + 1; |
394 |
> |
#else |
395 |
> |
exI = tempI + 1; |
396 |
> |
exJ = tempJ + 1; |
397 |
> |
exK = tempK + 1; |
398 |
> |
exL = tempL + 1; |
399 |
> |
#endif |
400 |
|
|
401 |
< |
info[k].excludes[j + excludeOffset]->setPair(exI, exJ); |
402 |
< |
#else // isn't MPI |
403 |
< |
info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1)); |
404 |
< |
#endif //is_mpi |
401 |
> |
info[k].excludes->addPair(exI, exJ); |
402 |
> |
info[k].excludes->addPair(exI, exK); |
403 |
> |
info[k].excludes->addPair(exI, exL); |
404 |
> |
info[k].excludes->addPair(exJ, exK); |
405 |
> |
info[k].excludes->addPair(exJ, exL); |
406 |
> |
info[k].excludes->addPair(exK, exL); |
407 |
|
} |
381 |
– |
excludeOffset += molInfo.nTorsions; |
408 |
|
|
409 |
+ |
for (j = 0; j < molInfo.nRigidBodies; j++){ |
410 |
|
|
411 |
< |
// send the arrays off to the forceField for init. |
411 |
> |
currentRigidBody = comp_stamps[stampID]->getRigidBody(j); |
412 |
> |
nMembers = currentRigidBody->getNMembers(); |
413 |
> |
|
414 |
> |
// Create the Rigid Body: |
415 |
> |
|
416 |
> |
myRB = new RigidBody(); |
417 |
> |
|
418 |
> |
for (rb1 = 0; rb1 < nMembers; rb1++) { |
419 |
> |
|
420 |
> |
// molI is atom numbering inside this molecule |
421 |
> |
molI = currentRigidBody->getMember(rb1); |
422 |
> |
|
423 |
> |
// tempI is atom numbering on local processor |
424 |
> |
tempI = molI + atomOffset; |
425 |
> |
|
426 |
> |
// currentAtom is the AtomStamp (which we need for |
427 |
> |
// rigid body reference positions) |
428 |
> |
currentAtom = comp_stamps[stampID]->getAtom(molI); |
429 |
|
|
430 |
+ |
// When we add to the rigid body, add the atom itself and |
431 |
+ |
// the stamp info: |
432 |
+ |
|
433 |
+ |
myRB->addAtom(info[k].atoms[tempI], currentAtom); |
434 |
+ |
|
435 |
+ |
// Add this atom to the Skip List for the integrators |
436 |
+ |
#ifdef IS_MPI |
437 |
+ |
slI = info[k].atoms[tempI]->getGlobalIndex(); |
438 |
+ |
#else |
439 |
+ |
slI = tempI; |
440 |
+ |
#endif |
441 |
+ |
skipList.insert(slI); |
442 |
+ |
|
443 |
+ |
} |
444 |
+ |
|
445 |
+ |
for(rb1 = 0; rb1 < nMembers - 1; rb1++) { |
446 |
+ |
for(rb2 = rb1+1; rb2 < nMembers; rb2++) { |
447 |
+ |
|
448 |
+ |
tempI = currentRigidBody->getMember(rb1); |
449 |
+ |
tempJ = currentRigidBody->getMember(rb2); |
450 |
+ |
|
451 |
+ |
// Some explanation is required here. |
452 |
+ |
// Fortran indexing starts at 1, while c indexing starts at 0 |
453 |
+ |
// Also, in parallel computations, the GlobalIndex is |
454 |
+ |
// used for the exclude list: |
455 |
+ |
|
456 |
+ |
#ifdef IS_MPI |
457 |
+ |
exI = info[k].atoms[tempI]->getGlobalIndex() + 1; |
458 |
+ |
exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1; |
459 |
+ |
#else |
460 |
+ |
exI = tempI + 1; |
461 |
+ |
exJ = tempJ + 1; |
462 |
+ |
#endif |
463 |
+ |
|
464 |
+ |
info[k].excludes->addPair(exI, exJ); |
465 |
+ |
|
466 |
+ |
} |
467 |
+ |
} |
468 |
+ |
} |
469 |
+ |
|
470 |
+ |
// send the arrays off to the forceField for init. |
471 |
+ |
|
472 |
|
the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms); |
473 |
|
the_ff->initializeBonds(molInfo.nBonds, molInfo.myBonds, theBonds); |
474 |
|
the_ff->initializeBends(molInfo.nBends, molInfo.myBends, theBends); |
475 |
|
the_ff->initializeTorsions(molInfo.nTorsions, molInfo.myTorsions, |
476 |
|
theTorsions); |
477 |
|
|
392 |
– |
|
478 |
|
info[k].molecules[i].initialize(molInfo); |
479 |
|
|
480 |
|
|
483 |
|
delete[] theBends; |
484 |
|
delete[] theTorsions; |
485 |
|
} |
486 |
+ |
|
487 |
+ |
// build up the integrableObjects vector: |
488 |
+ |
|
489 |
+ |
for (i = 0; i < info[k].n_atoms; i++) { |
490 |
+ |
|
491 |
+ |
#ifdef IS_MPI |
492 |
+ |
slI = info[k].atoms[i]->getGlobalIndex(); |
493 |
+ |
#else |
494 |
+ |
slI = i; |
495 |
+ |
#endif |
496 |
+ |
|
497 |
+ |
if (skipList.find(slI) == skipList.end()) { |
498 |
+ |
mySD = (StuntDouble *) info[k].atoms[i]; |
499 |
+ |
info[k].integrableObjects.push_back(mySD); |
500 |
+ |
} |
501 |
+ |
} |
502 |
+ |
for (i = 0; i < info[k].rigidBodies.size(); i++) { |
503 |
+ |
mySD = (StuntDouble *) info[k].rigidBodies[i]; |
504 |
+ |
info[k].integrableObjects.push_back(mySD); |
505 |
+ |
} |
506 |
+ |
|
507 |
|
} |
508 |
|
|
509 |
|
#ifdef IS_MPI |
513 |
|
|
514 |
|
// clean up the forcefield |
515 |
|
|
516 |
< |
the_ff->calcRcut(); |
516 |
> |
if (!globals->haveLJrcut()){ |
517 |
> |
|
518 |
> |
the_ff->calcRcut(); |
519 |
> |
|
520 |
> |
} else { |
521 |
> |
|
522 |
> |
the_ff->setRcut( globals->getLJrcut() ); |
523 |
> |
} |
524 |
> |
|
525 |
|
the_ff->cleanMe(); |
526 |
|
} |
527 |
|
|
722 |
|
} |
723 |
|
else if (!strcasecmp(force_field, "EAM")){ |
724 |
|
ffCase = FF_EAM; |
725 |
+ |
} |
726 |
+ |
else if (!strcasecmp(force_field, "WATER")){ |
727 |
+ |
ffCase = FF_H2O; |
728 |
|
} |
729 |
|
else{ |
730 |
|
sprintf(painCave.errMsg, "SimSetup Error. Unrecognized force field -> %s\n", |
928 |
|
for (int i = 0; i < nInfo; i++){ |
929 |
|
info[i].setSeed(seedValue); |
930 |
|
} |
931 |
< |
|
931 |
> |
|
932 |
|
#ifdef IS_MPI |
933 |
< |
strcpy(checkPointMsg, "Succesfully gathered all information from Bass\n"); |
933 |
> |
strcpy(checkPointMsg, "Successfully gathered all information from Bass\n"); |
934 |
|
MPIcheckPoint(); |
935 |
|
#endif // is_mpi |
936 |
|
} |
1227 |
|
the_ff = new EAM_FF(); |
1228 |
|
break; |
1229 |
|
|
1230 |
+ |
case FF_H2O: |
1231 |
+ |
the_ff = new WATER(); |
1232 |
+ |
break; |
1233 |
+ |
|
1234 |
|
default: |
1235 |
|
sprintf(painCave.errMsg, |
1236 |
|
"SimSetup Error. Unrecognized force field in case statement.\n"); |
1303 |
|
tot_bonds = 0; |
1304 |
|
tot_bends = 0; |
1305 |
|
tot_torsions = 0; |
1306 |
+ |
tot_rigid = 0; |
1307 |
|
for (i = 0; i < n_components; i++){ |
1308 |
|
tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms(); |
1309 |
|
tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds(); |
1310 |
|
tot_bends += components_nmol[i] * comp_stamps[i]->getNBends(); |
1311 |
|
tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions(); |
1312 |
+ |
tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies(); |
1313 |
|
} |
1314 |
< |
|
1314 |
> |
|
1315 |
|
tot_SRI = tot_bonds + tot_bends + tot_torsions; |
1316 |
|
molMembershipArray = new int[tot_atoms]; |
1317 |
|
|
1333 |
|
int i, j, k; |
1334 |
|
int localMol, allMol; |
1335 |
|
int local_atoms, local_bonds, local_bends, local_torsions, local_SRI; |
1336 |
+ |
int local_rigid; |
1337 |
|
|
1338 |
|
mpiSim = new mpiSimulation(info); |
1339 |
|
|
1350 |
|
local_bonds = 0; |
1351 |
|
local_bends = 0; |
1352 |
|
local_torsions = 0; |
1353 |
+ |
local_rigid = 0; |
1354 |
|
globalAtomIndex = 0; |
1355 |
|
|
1231 |
– |
|
1356 |
|
for (i = 0; i < n_components; i++){ |
1357 |
|
for (j = 0; j < components_nmol[i]; j++){ |
1358 |
|
if (mol2proc[allMol] == worldRank){ |
1360 |
|
local_bonds += comp_stamps[i]->getNBonds(); |
1361 |
|
local_bends += comp_stamps[i]->getNBends(); |
1362 |
|
local_torsions += comp_stamps[i]->getNTorsions(); |
1363 |
+ |
local_rigid += comp_stamps[i]->getNRigidBodies(); |
1364 |
|
localMol++; |
1365 |
|
} |
1366 |
|
for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){ |
1374 |
|
local_SRI = local_bonds + local_bends + local_torsions; |
1375 |
|
|
1376 |
|
info[0].n_atoms = mpiSim->getMyNlocal(); |
1377 |
+ |
|
1378 |
|
|
1379 |
|
if (local_atoms != info[0].n_atoms){ |
1380 |
|
sprintf(painCave.errMsg, |
1407 |
|
|
1408 |
|
Atom** the_atoms; |
1409 |
|
Molecule* the_molecules; |
1284 |
– |
Exclude** the_excludes; |
1410 |
|
|
1286 |
– |
|
1411 |
|
for (l = 0; l < nInfo; l++){ |
1412 |
|
// create the atom and short range interaction arrays |
1413 |
|
|
1450 |
|
|
1451 |
|
#endif // is_mpi |
1452 |
|
|
1453 |
< |
|
1454 |
< |
if (info[l].n_SRI){ |
1455 |
< |
Exclude::createArray(info[l].n_SRI); |
1332 |
< |
the_excludes = new Exclude * [info[l].n_SRI]; |
1333 |
< |
for (int ex = 0; ex < info[l].n_SRI; ex++){ |
1334 |
< |
the_excludes[ex] = new Exclude(ex); |
1335 |
< |
} |
1336 |
< |
info[l].globalExcludes = new int; |
1337 |
< |
info[l].n_exclude = info[l].n_SRI; |
1338 |
< |
} |
1339 |
< |
else{ |
1340 |
< |
Exclude::createArray(1); |
1341 |
< |
the_excludes = new Exclude * ; |
1342 |
< |
the_excludes[0] = new Exclude(0); |
1343 |
< |
the_excludes[0]->setPair(0, 0); |
1344 |
< |
info[l].globalExcludes = new int; |
1345 |
< |
info[l].globalExcludes[0] = 0; |
1346 |
< |
info[l].n_exclude = 0; |
1347 |
< |
} |
1348 |
< |
|
1453 |
> |
info[l].globalExcludes = new int; |
1454 |
> |
info[l].globalExcludes[0] = 0; |
1455 |
> |
|
1456 |
|
// set the arrays into the SimInfo object |
1457 |
|
|
1458 |
|
info[l].atoms = the_atoms; |
1459 |
|
info[l].molecules = the_molecules; |
1460 |
|
info[l].nGlobalExcludes = 0; |
1354 |
– |
info[l].excludes = the_excludes; |
1461 |
|
|
1462 |
|
the_ff->setSimInfo(info); |
1463 |
|
} |
1731 |
|
|
1732 |
|
theInfo.addProperty(zconsForcePolicy); |
1733 |
|
|
1734 |
+ |
//set zcons gap |
1735 |
+ |
DoubleData* zconsGap = new DoubleData(); |
1736 |
+ |
zconsGap->setID(ZCONSGAP_ID); |
1737 |
+ |
|
1738 |
+ |
if (globals->haveZConsGap()){ |
1739 |
+ |
zconsGap->setData(globals->getZconsGap()); |
1740 |
+ |
theInfo.addProperty(zconsGap); |
1741 |
+ |
} |
1742 |
+ |
|
1743 |
+ |
//set zcons fixtime |
1744 |
+ |
DoubleData* zconsFixtime = new DoubleData(); |
1745 |
+ |
zconsFixtime->setID(ZCONSFIXTIME_ID); |
1746 |
+ |
|
1747 |
+ |
if (globals->haveZConsFixTime()){ |
1748 |
+ |
zconsFixtime->setData(globals->getZconsFixtime()); |
1749 |
+ |
theInfo.addProperty(zconsFixtime); |
1750 |
+ |
} |
1751 |
+ |
|
1752 |
+ |
//set zconsUsingSMD |
1753 |
+ |
IntData* zconsUsingSMD = new IntData(); |
1754 |
+ |
zconsUsingSMD->setID(ZCONSUSINGSMD_ID); |
1755 |
+ |
|
1756 |
+ |
if (globals->haveZConsUsingSMD()){ |
1757 |
+ |
zconsUsingSMD->setData(globals->getZconsUsingSMD()); |
1758 |
+ |
theInfo.addProperty(zconsUsingSMD); |
1759 |
+ |
} |
1760 |
+ |
|
1761 |
|
//Determine the name of ouput file and add it into SimInfo's property list |
1762 |
|
//Be careful, do not use inFileName, since it is a pointer which |
1763 |
|
//point to a string at master node, and slave nodes do not contain that string |
1787 |
|
tempParaItem.zPos = zconStamp[i]->getZpos(); |
1788 |
|
tempParaItem.zconsIndex = zconStamp[i]->getMolIndex(); |
1789 |
|
tempParaItem.kRatio = zconStamp[i]->getKratio(); |
1790 |
< |
|
1790 |
> |
tempParaItem.havingCantVel = zconStamp[i]->haveCantVel(); |
1791 |
> |
tempParaItem.cantVel = zconStamp[i]->getCantVel(); |
1792 |
|
zconsParaData->addItem(tempParaItem); |
1793 |
|
} |
1794 |
|
|
1806 |
|
//push data into siminfo, therefore, we can retrieve later |
1807 |
|
theInfo.addProperty(zconsParaData); |
1808 |
|
} |
1809 |
+ |
|
1810 |
+ |
void SimSetup::makeMinimizer(){ |
1811 |
+ |
|
1812 |
+ |
OOPSEMinimizer* myOOPSEMinimizer; |
1813 |
+ |
MinimizerParameterSet* param; |
1814 |
+ |
char minimizerName[100]; |
1815 |
+ |
|
1816 |
+ |
for (int i = 0; i < nInfo; i++){ |
1817 |
+ |
|
1818 |
+ |
//prepare parameter set for minimizer |
1819 |
+ |
param = new MinimizerParameterSet(); |
1820 |
+ |
param->setDefaultParameter(); |
1821 |
+ |
|
1822 |
+ |
if (globals->haveMinimizer()){ |
1823 |
+ |
param->setFTol(globals->getMinFTol()); |
1824 |
+ |
} |
1825 |
+ |
|
1826 |
+ |
if (globals->haveMinGTol()){ |
1827 |
+ |
param->setGTol(globals->getMinGTol()); |
1828 |
+ |
} |
1829 |
+ |
|
1830 |
+ |
if (globals->haveMinMaxIter()){ |
1831 |
+ |
param->setMaxIteration(globals->getMinMaxIter()); |
1832 |
+ |
} |
1833 |
+ |
|
1834 |
+ |
if (globals->haveMinWriteFrq()){ |
1835 |
+ |
param->setMaxIteration(globals->getMinMaxIter()); |
1836 |
+ |
} |
1837 |
+ |
|
1838 |
+ |
if (globals->haveMinWriteFrq()){ |
1839 |
+ |
param->setWriteFrq(globals->getMinWriteFrq()); |
1840 |
+ |
} |
1841 |
+ |
|
1842 |
+ |
if (globals->haveMinStepSize()){ |
1843 |
+ |
param->setStepSize(globals->getMinStepSize()); |
1844 |
+ |
} |
1845 |
+ |
|
1846 |
+ |
if (globals->haveMinLSMaxIter()){ |
1847 |
+ |
param->setLineSearchMaxIteration(globals->getMinLSMaxIter()); |
1848 |
+ |
} |
1849 |
+ |
|
1850 |
+ |
if (globals->haveMinLSTol()){ |
1851 |
+ |
param->setLineSearchTol(globals->getMinLSTol()); |
1852 |
+ |
} |
1853 |
+ |
|
1854 |
+ |
strcpy(minimizerName, globals->getMinimizer()); |
1855 |
+ |
|
1856 |
+ |
if (!strcasecmp(minimizerName, "CG")){ |
1857 |
+ |
myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param); |
1858 |
+ |
} |
1859 |
+ |
else if (!strcasecmp(minimizerName, "SD")){ |
1860 |
+ |
//myOOPSEMinimizer = MinimizerFactory.creatMinimizer("", &(info[i]), the_ff, param); |
1861 |
+ |
myOOPSEMinimizer = new SDMinimizer(&(info[i]), the_ff, param); |
1862 |
+ |
} |
1863 |
+ |
else{ |
1864 |
+ |
sprintf(painCave.errMsg, |
1865 |
+ |
"SimSetup error: Unrecognized Minimizer, use Conjugate Gradient \n"); |
1866 |
+ |
painCave.isFatal = 0; |
1867 |
+ |
simError(); |
1868 |
+ |
|
1869 |
+ |
myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param); |
1870 |
+ |
} |
1871 |
+ |
info[i].the_integrator = myOOPSEMinimizer; |
1872 |
+ |
|
1873 |
+ |
//store the minimizer into simInfo |
1874 |
+ |
info[i].the_minimizer = myOOPSEMinimizer; |
1875 |
+ |
info[i].has_minimizer = true; |
1876 |
+ |
} |
1877 |
+ |
|
1878 |
+ |
} |