1 |
|
#include <algorithm> |
2 |
< |
#include <cstdlib> |
2 |
> |
#include <stdlib.h> |
3 |
|
#include <iostream> |
4 |
< |
#include <cmath> |
4 |
> |
#include <math.h> |
5 |
|
#include <string> |
6 |
|
#include <sprng.h> |
7 |
– |
|
7 |
|
#include "SimSetup.hpp" |
8 |
|
#include "ReadWrite.hpp" |
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" |
24 |
|
#define NVT_ENS 1 |
25 |
|
#define NPTi_ENS 2 |
26 |
|
#define NPTf_ENS 3 |
27 |
< |
#define NPTim_ENS 4 |
26 |
< |
#define NPTfm_ENS 5 |
27 |
> |
#define NPTxyz_ENS 4 |
28 |
|
|
28 |
– |
#define FF_DUFF 0 |
29 |
– |
#define FF_LJ 1 |
30 |
– |
#define FF_EAM 2 |
29 |
|
|
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 |
|
|
37 |
+ |
/** |
38 |
+ |
* Check whether dividend is divisble by divisor or not |
39 |
+ |
*/ |
40 |
+ |
bool isDivisible(double dividend, double divisor){ |
41 |
+ |
double tolerance = 0.000001; |
42 |
+ |
double quotient; |
43 |
+ |
double diff; |
44 |
+ |
int intQuotient; |
45 |
+ |
|
46 |
+ |
quotient = dividend / divisor; |
47 |
+ |
|
48 |
+ |
if (quotient < 0) |
49 |
+ |
quotient = -quotient; |
50 |
+ |
|
51 |
+ |
intQuotient = int (quotient + tolerance); |
52 |
+ |
|
53 |
+ |
diff = fabs(fabs(dividend) - intQuotient * fabs(divisor)); |
54 |
+ |
|
55 |
+ |
if (diff <= tolerance) |
56 |
+ |
return true; |
57 |
+ |
else |
58 |
+ |
return false; |
59 |
+ |
} |
60 |
+ |
|
61 |
|
SimSetup::SimSetup(){ |
62 |
+ |
|
63 |
+ |
initSuspend = false; |
64 |
|
isInfoArray = 0; |
65 |
|
nInfo = 1; |
66 |
|
|
83 |
|
info = the_info; |
84 |
|
nInfo = theNinfo; |
85 |
|
isInfoArray = 1; |
86 |
+ |
initSuspend = true; |
87 |
|
} |
88 |
|
|
89 |
|
|
122 |
|
#endif // is_mpi |
123 |
|
|
124 |
|
void SimSetup::createSim(void){ |
95 |
– |
int i, j, k, globalAtomIndex; |
125 |
|
|
126 |
|
// gather all of the information from the Bass file |
127 |
|
|
137 |
|
|
138 |
|
// initialize the system coordinates |
139 |
|
|
140 |
< |
if (!isInfoArray){ |
140 |
> |
if ( !initSuspend ){ |
141 |
|
initSystemCoords(); |
142 |
+ |
|
143 |
+ |
if( !(globals->getUseInitTime()) ) |
144 |
+ |
info[0].currentTime = 0.0; |
145 |
|
} |
146 |
|
|
147 |
|
// make the output filenames |
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, l; |
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 |
|
|
154 |
– |
|
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]); |
174 |
– |
molInfo.myExcludes = &(info[k].excludes[excludeOffset]); |
175 |
– |
molInfo.myBonds = new Bond * [molInfo.nBonds]; |
176 |
– |
molInfo.myBends = new Bend * [molInfo.nBends]; |
177 |
– |
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; |
201 |
< |
uy = uy / u; |
202 |
< |
uz = uz / u; |
203 |
< |
|
204 |
< |
dAtom->setSUx(ux); |
205 |
< |
dAtom->setSUy(uy); |
206 |
< |
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 |
|
|
230 |
– |
// exclude_I must always be the smaller of the pair |
231 |
– |
if (exI > exJ){ |
232 |
– |
tempEx = exI; |
233 |
– |
exI = exJ; |
234 |
– |
exJ = tempEx; |
235 |
– |
} |
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); |
243 |
< |
#else // isn't MPI |
244 |
< |
|
245 |
< |
info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1)); |
246 |
< |
#endif //is_mpi |
290 |
> |
info[k].excludes->addPair(exI, exJ); |
291 |
|
} |
248 |
– |
excludeOffset += molInfo.nBonds; |
292 |
|
|
293 |
|
//make the bends |
294 |
|
for (j = 0; j < molInfo.nBends; j++){ |
338 |
|
} |
339 |
|
} |
340 |
|
|
341 |
< |
if (!theBends[j].isGhost){ |
342 |
< |
exI = theBends[j].a; |
343 |
< |
exJ = theBends[j].c; |
344 |
< |
} |
345 |
< |
else{ |
303 |
< |
exI = theBends[j].a; |
304 |
< |
exJ = theBends[j].b; |
305 |
< |
} |
306 |
< |
|
307 |
< |
// exclude_I must always be the smaller of the pair |
308 |
< |
if (exI > exJ){ |
309 |
< |
tempEx = exI; |
310 |
< |
exI = exJ; |
311 |
< |
exJ = tempEx; |
312 |
< |
} |
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 |
|
} |
324 |
– |
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 |
|
|
336 |
– |
// exclude_I must always be the smaller of the pair |
337 |
– |
if (exI > exJ){ |
338 |
– |
tempEx = exI; |
339 |
– |
exI = exJ; |
340 |
– |
exJ = tempEx; |
341 |
– |
} |
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 |
|
} |
353 |
– |
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 |
|
|
364 |
– |
|
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 |
|
|
695 |
|
|
696 |
|
|
697 |
|
void SimSetup::gatherInfo(void){ |
698 |
< |
int i, j, k; |
698 |
> |
int i; |
699 |
|
|
700 |
|
ensembleCase = -1; |
701 |
|
ffCase = -1; |
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", |
731 |
|
force_field); |
748 |
|
} |
749 |
|
else if (!strcasecmp(ensemble, "NPTf")){ |
750 |
|
ensembleCase = NPTf_ENS; |
606 |
– |
} |
607 |
– |
else if (!strcasecmp(ensemble, "NPTim")){ |
608 |
– |
ensembleCase = NPTim_ENS; |
751 |
|
} |
752 |
< |
else if (!strcasecmp(ensemble, "NPTfm")){ |
753 |
< |
ensembleCase = NPTfm_ENS; |
752 |
> |
else if (!strcasecmp(ensemble, "NPTxyz")){ |
753 |
> |
ensembleCase = NPTxyz_ENS; |
754 |
|
} |
755 |
|
else{ |
756 |
|
sprintf(painCave.errMsg, |
757 |
< |
"SimSetup Warning. Unrecognized Ensemble -> %s, " |
758 |
< |
"reverting to NVE for this simulation.\n", |
757 |
> |
"SimSetup Warning. Unrecognized Ensemble -> %s \n" |
758 |
> |
"\treverting to NVE for this simulation.\n", |
759 |
|
ensemble); |
760 |
|
painCave.isFatal = 0; |
761 |
|
simError(); |
787 |
|
if (!the_components[i]->haveNMol()){ |
788 |
|
// we have a problem |
789 |
|
sprintf(painCave.errMsg, |
790 |
< |
"SimSetup Error. No global NMol or component NMol" |
791 |
< |
" given. Cannot calculate the number of atoms.\n"); |
790 |
> |
"SimSetup Error. No global NMol or component NMol given.\n" |
791 |
> |
"\tCannot calculate the number of atoms.\n"); |
792 |
|
painCave.isFatal = 1; |
793 |
|
simError(); |
794 |
|
} |
808 |
|
simError(); |
809 |
|
} |
810 |
|
|
811 |
+ |
//check whether sample time, status time, thermal time and reset time are divisble by dt |
812 |
+ |
if (!isDivisible(globals->getSampleTime(), globals->getDt())){ |
813 |
+ |
sprintf(painCave.errMsg, |
814 |
+ |
"Sample time is not divisible by dt.\n" |
815 |
+ |
"\tThis will result in samples that are not uniformly\n" |
816 |
+ |
"\tdistributed in time. If this is a problem, change\n" |
817 |
+ |
"\tyour sampleTime variable.\n"); |
818 |
+ |
painCave.isFatal = 0; |
819 |
+ |
simError(); |
820 |
+ |
} |
821 |
+ |
|
822 |
+ |
if (globals->haveStatusTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){ |
823 |
+ |
sprintf(painCave.errMsg, |
824 |
+ |
"Status time is not divisible by dt.\n" |
825 |
+ |
"\tThis will result in status reports that are not uniformly\n" |
826 |
+ |
"\tdistributed in time. If this is a problem, change \n" |
827 |
+ |
"\tyour statusTime variable.\n"); |
828 |
+ |
painCave.isFatal = 0; |
829 |
+ |
simError(); |
830 |
+ |
} |
831 |
+ |
|
832 |
+ |
if (globals->haveThermalTime() && !isDivisible(globals->getThermalTime(), globals->getDt())){ |
833 |
+ |
sprintf(painCave.errMsg, |
834 |
+ |
"Thermal time is not divisible by dt.\n" |
835 |
+ |
"\tThis will result in thermalizations that are not uniformly\n" |
836 |
+ |
"\tdistributed in time. If this is a problem, change \n" |
837 |
+ |
"\tyour thermalTime variable.\n"); |
838 |
+ |
painCave.isFatal = 0; |
839 |
+ |
simError(); |
840 |
+ |
} |
841 |
+ |
|
842 |
+ |
if (globals->haveResetTime() && !isDivisible(globals->getResetTime(), globals->getDt())){ |
843 |
+ |
sprintf(painCave.errMsg, |
844 |
+ |
"Reset time is not divisible by dt.\n" |
845 |
+ |
"\tThis will result in integrator resets that are not uniformly\n" |
846 |
+ |
"\tdistributed in time. If this is a problem, change\n" |
847 |
+ |
"\tyour resetTime variable.\n"); |
848 |
+ |
painCave.isFatal = 0; |
849 |
+ |
simError(); |
850 |
+ |
} |
851 |
+ |
|
852 |
|
// set the status, sample, and thermal kick times |
853 |
|
|
854 |
|
for (i = 0; i < nInfo; i++){ |
878 |
|
} |
879 |
|
|
880 |
|
// check for the temperature set flag |
881 |
< |
|
881 |
> |
|
882 |
|
if (globals->haveTempSet()) |
883 |
|
info[i].setTemp = globals->getTempSet(); |
884 |
|
|
885 |
< |
// get some of the tricky things that may still be in the globals |
703 |
< |
|
704 |
< |
double boxVector[3]; |
705 |
< |
if (globals->haveBox()){ |
706 |
< |
boxVector[0] = globals->getBox(); |
707 |
< |
boxVector[1] = globals->getBox(); |
708 |
< |
boxVector[2] = globals->getBox(); |
885 |
> |
// check for the extended State init |
886 |
|
|
887 |
< |
info[i].setBox(boxVector); |
888 |
< |
} |
889 |
< |
else if (globals->haveDensity()){ |
713 |
< |
double vol; |
714 |
< |
vol = (double) tot_nmol / globals->getDensity(); |
715 |
< |
boxVector[0] = pow(vol, (1.0 / 3.0)); |
716 |
< |
boxVector[1] = boxVector[0]; |
717 |
< |
boxVector[2] = boxVector[0]; |
718 |
< |
|
719 |
< |
info[i].setBox(boxVector); |
720 |
< |
} |
721 |
< |
else{ |
722 |
< |
if (!globals->haveBoxX()){ |
723 |
< |
sprintf(painCave.errMsg, |
724 |
< |
"SimSetup error, no periodic BoxX size given.\n"); |
725 |
< |
painCave.isFatal = 1; |
726 |
< |
simError(); |
727 |
< |
} |
728 |
< |
boxVector[0] = globals->getBoxX(); |
729 |
< |
|
730 |
< |
if (!globals->haveBoxY()){ |
731 |
< |
sprintf(painCave.errMsg, |
732 |
< |
"SimSetup error, no periodic BoxY size given.\n"); |
733 |
< |
painCave.isFatal = 1; |
734 |
< |
simError(); |
735 |
< |
} |
736 |
< |
boxVector[1] = globals->getBoxY(); |
737 |
< |
|
738 |
< |
if (!globals->haveBoxZ()){ |
739 |
< |
sprintf(painCave.errMsg, |
740 |
< |
"SimSetup error, no periodic BoxZ size given.\n"); |
741 |
< |
painCave.isFatal = 1; |
742 |
< |
simError(); |
743 |
< |
} |
744 |
< |
boxVector[2] = globals->getBoxZ(); |
745 |
< |
|
746 |
< |
info[i].setBox(boxVector); |
747 |
< |
} |
887 |
> |
info[i].useInitXSstate = globals->getUseInitXSstate(); |
888 |
> |
info[i].orthoTolerance = globals->getOrthoBoxTolerance(); |
889 |
> |
|
890 |
|
} |
891 |
< |
|
891 |
> |
|
892 |
|
//setup seed for random number generator |
893 |
|
int seedValue; |
894 |
|
|
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 |
|
} |
963 |
|
|
964 |
|
if (!globals->haveECR()){ |
965 |
|
sprintf(painCave.errMsg, |
966 |
< |
"SimSetup Warning: using default value of 1/2 the smallest " |
967 |
< |
"box length for the electrostaticCutoffRadius.\n" |
968 |
< |
"I hope you have a very fast processor!\n"); |
966 |
> |
"SimSetup Warning: No value was set for electrostaticCutoffRadius.\n" |
967 |
> |
"\tOOPSE will use a default value of 15.0 angstroms" |
968 |
> |
"\tfor the electrostaticCutoffRadius.\n"); |
969 |
|
painCave.isFatal = 0; |
970 |
|
simError(); |
971 |
< |
double smallest; |
830 |
< |
smallest = info[i].boxL[0]; |
831 |
< |
if (info[i].boxL[1] <= smallest) |
832 |
< |
smallest = info[i].boxL[1]; |
833 |
< |
if (info[i].boxL[2] <= smallest) |
834 |
< |
smallest = info[i].boxL[2]; |
835 |
< |
theEcr = 0.5 * smallest; |
971 |
> |
theEcr = 15.0; |
972 |
|
} |
973 |
|
else{ |
974 |
|
theEcr = globals->getECR(); |
976 |
|
|
977 |
|
if (!globals->haveEST()){ |
978 |
|
sprintf(painCave.errMsg, |
979 |
< |
"SimSetup Warning: using default value of 0.05 * the " |
980 |
< |
"electrostaticCutoffRadius for the electrostaticSkinThickness\n"); |
979 |
> |
"SimSetup Warning: No value was set for electrostaticSkinThickness.\n" |
980 |
> |
"\tOOPSE will use a default value of\n" |
981 |
> |
"\t0.05 * electrostaticCutoffRadius\n" |
982 |
> |
"\tfor the electrostaticSkinThickness\n"); |
983 |
|
painCave.isFatal = 0; |
984 |
|
simError(); |
985 |
|
theEst = 0.05 * theEcr; |
988 |
|
theEst = globals->getEST(); |
989 |
|
} |
990 |
|
|
991 |
< |
info[i].setEcr(theEcr, theEst); |
991 |
> |
info[i].setDefaultEcr(theEcr, theEst); |
992 |
|
|
993 |
|
if (!globals->haveDielectric()){ |
994 |
|
sprintf(painCave.errMsg, |
995 |
< |
"SimSetup Error: You are trying to use Reaction Field without" |
996 |
< |
"setting a dielectric constant!\n"); |
995 |
> |
"SimSetup Error: No Dielectric constant was set.\n" |
996 |
> |
"\tYou are trying to use Reaction Field without" |
997 |
> |
"\tsetting a dielectric constant!\n"); |
998 |
|
painCave.isFatal = 1; |
999 |
|
simError(); |
1000 |
|
} |
1004 |
|
if (usesDipoles){ |
1005 |
|
if (!globals->haveECR()){ |
1006 |
|
sprintf(painCave.errMsg, |
1007 |
< |
"SimSetup Warning: using default value of 1/2 the smallest " |
1008 |
< |
"box length for the electrostaticCutoffRadius.\n" |
1009 |
< |
"I hope you have a very fast processor!\n"); |
1010 |
< |
painCave.isFatal = 0; |
1011 |
< |
simError(); |
1012 |
< |
double smallest; |
874 |
< |
smallest = info[i].boxL[0]; |
875 |
< |
if (info[i].boxL[1] <= smallest) |
876 |
< |
smallest = info[i].boxL[1]; |
877 |
< |
if (info[i].boxL[2] <= smallest) |
878 |
< |
smallest = info[i].boxL[2]; |
879 |
< |
theEcr = 0.5 * smallest; |
1007 |
> |
"SimSetup Warning: No value was set for electrostaticCutoffRadius.\n" |
1008 |
> |
"\tOOPSE will use a default value of 15.0 angstroms" |
1009 |
> |
"\tfor the electrostaticCutoffRadius.\n"); |
1010 |
> |
painCave.isFatal = 0; |
1011 |
> |
simError(); |
1012 |
> |
theEcr = 15.0; |
1013 |
|
} |
1014 |
|
else{ |
1015 |
|
theEcr = globals->getECR(); |
1016 |
|
} |
1017 |
< |
|
1017 |
> |
|
1018 |
|
if (!globals->haveEST()){ |
1019 |
|
sprintf(painCave.errMsg, |
1020 |
< |
"SimSetup Warning: using default value of 0.05 * the " |
1021 |
< |
"electrostaticCutoffRadius for the " |
1022 |
< |
"electrostaticSkinThickness\n"); |
1020 |
> |
"SimSetup Warning: No value was set for electrostaticSkinThickness.\n" |
1021 |
> |
"\tOOPSE will use a default value of\n" |
1022 |
> |
"\t0.05 * electrostaticCutoffRadius\n" |
1023 |
> |
"\tfor the electrostaticSkinThickness\n"); |
1024 |
|
painCave.isFatal = 0; |
1025 |
|
simError(); |
1026 |
|
theEst = 0.05 * theEcr; |
1028 |
|
else{ |
1029 |
|
theEst = globals->getEST(); |
1030 |
|
} |
1031 |
< |
|
1032 |
< |
info[i].setEcr(theEcr, theEst); |
1031 |
> |
|
1032 |
> |
info[i].setDefaultEcr(theEcr, theEst); |
1033 |
|
} |
1034 |
|
} |
1035 |
|
} |
902 |
– |
|
1036 |
|
#ifdef IS_MPI |
1037 |
|
strcpy(checkPointMsg, "post processing checks out"); |
1038 |
|
MPIcheckPoint(); |
1039 |
|
#endif // is_mpi |
1040 |
|
} |
1041 |
< |
|
1041 |
> |
|
1042 |
|
void SimSetup::initSystemCoords(void){ |
1043 |
|
int i; |
1044 |
|
|
1055 |
|
if (worldRank == 0){ |
1056 |
|
#endif //is_mpi |
1057 |
|
inName = globals->getInitialConfig(); |
925 |
– |
double* tempDouble = new double[1000000]; |
1058 |
|
fileInit = new InitializeFromFile(inName); |
1059 |
|
#ifdef IS_MPI |
1060 |
|
} |
1066 |
|
delete fileInit; |
1067 |
|
} |
1068 |
|
else{ |
1069 |
< |
#ifdef IS_MPI |
938 |
< |
|
1069 |
> |
|
1070 |
|
// no init from bass |
1071 |
< |
|
1071 |
> |
|
1072 |
|
sprintf(painCave.errMsg, |
1073 |
< |
"Cannot intialize a parallel simulation without an initial configuration file.\n"); |
1074 |
< |
painCave.isFatal; |
1073 |
> |
"Cannot intialize a simulation without an initial configuration file.\n"); |
1074 |
> |
painCave.isFatal = 1;; |
1075 |
|
simError(); |
1076 |
< |
|
946 |
< |
#else |
947 |
< |
|
948 |
< |
initFromBass(); |
949 |
< |
|
950 |
< |
|
951 |
< |
#endif |
1076 |
> |
|
1077 |
|
} |
1078 |
|
|
1079 |
|
#ifdef IS_MPI |
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"); |
1295 |
|
} |
1296 |
|
|
1297 |
|
void SimSetup::calcSysValues(void){ |
1298 |
< |
int i, j, k; |
1298 |
> |
int i; |
1299 |
|
|
1300 |
|
int* molMembershipArray; |
1301 |
|
|
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 |
|
|
1223 |
– |
|
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, |
1381 |
< |
"SimSetup error: mpiSim's localAtom (%d) and SimSetup's" |
1382 |
< |
" localAtom (%d) are not equal.\n", |
1381 |
> |
"SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n" |
1382 |
> |
"\tlocalAtom (%d) are not equal.\n", |
1383 |
|
info[0].n_atoms, local_atoms); |
1384 |
|
painCave.isFatal = 1; |
1385 |
|
simError(); |
1399 |
|
|
1400 |
|
|
1401 |
|
void SimSetup::makeSysArrays(void){ |
1402 |
< |
int i, j, k, l; |
1402 |
> |
|
1403 |
> |
#ifndef IS_MPI |
1404 |
> |
int k, j; |
1405 |
> |
#endif // is_mpi |
1406 |
> |
int i, l; |
1407 |
|
|
1408 |
|
Atom** the_atoms; |
1409 |
|
Molecule* the_molecules; |
1272 |
– |
Exclude** the_excludes; |
1410 |
|
|
1274 |
– |
|
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); |
1320 |
< |
the_excludes = new Exclude * [info[l].n_SRI]; |
1321 |
< |
for (int ex = 0; ex < info[l].n_SRI; ex++){ |
1322 |
< |
the_excludes[ex] = new Exclude(ex); |
1323 |
< |
} |
1324 |
< |
info[l].globalExcludes = new int; |
1325 |
< |
info[l].n_exclude = info[l].n_SRI; |
1326 |
< |
} |
1327 |
< |
else{ |
1328 |
< |
Exclude::createArray(1); |
1329 |
< |
the_excludes = new Exclude * ; |
1330 |
< |
the_excludes[0] = new Exclude(0); |
1331 |
< |
the_excludes[0]->setPair(0, 0); |
1332 |
< |
info[l].globalExcludes = new int; |
1333 |
< |
info[l].globalExcludes[0] = 0; |
1334 |
< |
info[l].n_exclude = 0; |
1335 |
< |
} |
1336 |
< |
|
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; |
1342 |
– |
info[l].excludes = the_excludes; |
1461 |
|
|
1462 |
|
the_ff->setSimInfo(info); |
1463 |
|
} |
1466 |
|
void SimSetup::makeIntegrator(void){ |
1467 |
|
int k; |
1468 |
|
|
1469 |
+ |
NVE<RealIntegrator>* myNVE = NULL; |
1470 |
|
NVT<RealIntegrator>* myNVT = NULL; |
1471 |
< |
NPTi<RealIntegrator>* myNPTi = NULL; |
1472 |
< |
NPTf<RealIntegrator>* myNPTf = NULL; |
1473 |
< |
NPTim<RealIntegrator>* myNPTim = NULL; |
1355 |
< |
NPTfm<RealIntegrator>* myNPTfm = NULL; |
1471 |
> |
NPTi<NPT<RealIntegrator> >* myNPTi = NULL; |
1472 |
> |
NPTf<NPT<RealIntegrator> >* myNPTf = NULL; |
1473 |
> |
NPTxyz<NPT<RealIntegrator> >* myNPTxyz = NULL; |
1474 |
|
|
1475 |
|
for (k = 0; k < nInfo; k++){ |
1476 |
|
switch (ensembleCase){ |
1477 |
|
case NVE_ENS: |
1478 |
|
if (globals->haveZconstraints()){ |
1479 |
|
setupZConstraint(info[k]); |
1480 |
< |
new ZConstraint<NVE<RealIntegrator> >(&(info[k]), the_ff); |
1480 |
> |
myNVE = new ZConstraint<NVE<RealIntegrator> >(&(info[k]), the_ff); |
1481 |
|
} |
1482 |
< |
else |
1483 |
< |
new NVE<RealIntegrator>(&(info[k]), the_ff); |
1482 |
> |
else{ |
1483 |
> |
myNVE = new NVE<RealIntegrator>(&(info[k]), the_ff); |
1484 |
> |
} |
1485 |
> |
|
1486 |
> |
info->the_integrator = myNVE; |
1487 |
|
break; |
1488 |
|
|
1489 |
|
case NVT_ENS: |
1501 |
|
else{ |
1502 |
|
sprintf(painCave.errMsg, |
1503 |
|
"SimSetup error: If you use the NVT\n" |
1504 |
< |
" ensemble, you must set tauThermostat.\n"); |
1504 |
> |
"\tensemble, you must set tauThermostat.\n"); |
1505 |
|
painCave.isFatal = 1; |
1506 |
|
simError(); |
1507 |
|
} |
1508 |
+ |
|
1509 |
+ |
info->the_integrator = myNVT; |
1510 |
|
break; |
1511 |
|
|
1512 |
|
case NPTi_ENS: |
1513 |
|
if (globals->haveZconstraints()){ |
1514 |
|
setupZConstraint(info[k]); |
1515 |
< |
myNPTi = new ZConstraint<NPTi<RealIntegrator> >(&(info[k]), the_ff); |
1515 |
> |
myNPTi = new ZConstraint<NPTi<NPT <RealIntegrator> > >(&(info[k]), the_ff); |
1516 |
|
} |
1517 |
|
else |
1518 |
< |
myNPTi = new NPTi<RealIntegrator>(&(info[k]), the_ff); |
1518 |
> |
myNPTi = new NPTi<NPT<RealIntegrator> >(&(info[k]), the_ff); |
1519 |
|
|
1520 |
|
myNPTi->setTargetTemp(globals->getTargetTemp()); |
1521 |
|
|
1524 |
|
else{ |
1525 |
|
sprintf(painCave.errMsg, |
1526 |
|
"SimSetup error: If you use a constant pressure\n" |
1527 |
< |
" ensemble, you must set targetPressure in the BASS file.\n"); |
1527 |
> |
"\tensemble, you must set targetPressure in the BASS file.\n"); |
1528 |
|
painCave.isFatal = 1; |
1529 |
|
simError(); |
1530 |
|
} |
1534 |
|
else{ |
1535 |
|
sprintf(painCave.errMsg, |
1536 |
|
"SimSetup error: If you use an NPT\n" |
1537 |
< |
" ensemble, you must set tauThermostat.\n"); |
1537 |
> |
"\tensemble, you must set tauThermostat.\n"); |
1538 |
|
painCave.isFatal = 1; |
1539 |
|
simError(); |
1540 |
|
} |
1544 |
|
else{ |
1545 |
|
sprintf(painCave.errMsg, |
1546 |
|
"SimSetup error: If you use an NPT\n" |
1547 |
< |
" ensemble, you must set tauBarostat.\n"); |
1547 |
> |
"\tensemble, you must set tauBarostat.\n"); |
1548 |
|
painCave.isFatal = 1; |
1549 |
|
simError(); |
1550 |
|
} |
1551 |
+ |
|
1552 |
+ |
info->the_integrator = myNPTi; |
1553 |
|
break; |
1554 |
|
|
1555 |
|
case NPTf_ENS: |
1556 |
|
if (globals->haveZconstraints()){ |
1557 |
|
setupZConstraint(info[k]); |
1558 |
< |
myNPTf = new ZConstraint<NPTf<RealIntegrator> >(&(info[k]), the_ff); |
1558 |
> |
myNPTf = new ZConstraint<NPTf<NPT <RealIntegrator> > >(&(info[k]), the_ff); |
1559 |
|
} |
1560 |
|
else |
1561 |
< |
myNPTf = new NPTf<RealIntegrator>(&(info[k]), the_ff); |
1561 |
> |
myNPTf = new NPTf<NPT <RealIntegrator> >(&(info[k]), the_ff); |
1562 |
|
|
1563 |
|
myNPTf->setTargetTemp(globals->getTargetTemp()); |
1564 |
|
|
1567 |
|
else{ |
1568 |
|
sprintf(painCave.errMsg, |
1569 |
|
"SimSetup error: If you use a constant pressure\n" |
1570 |
< |
" ensemble, you must set targetPressure in the BASS file.\n"); |
1570 |
> |
"\tensemble, you must set targetPressure in the BASS file.\n"); |
1571 |
|
painCave.isFatal = 1; |
1572 |
|
simError(); |
1573 |
|
} |
1574 |
|
|
1575 |
|
if (globals->haveTauThermostat()) |
1576 |
|
myNPTf->setTauThermostat(globals->getTauThermostat()); |
1577 |
+ |
|
1578 |
|
else{ |
1579 |
|
sprintf(painCave.errMsg, |
1580 |
|
"SimSetup error: If you use an NPT\n" |
1581 |
< |
" ensemble, you must set tauThermostat.\n"); |
1581 |
> |
"\tensemble, you must set tauThermostat.\n"); |
1582 |
|
painCave.isFatal = 1; |
1583 |
|
simError(); |
1584 |
|
} |
1585 |
|
|
1586 |
|
if (globals->haveTauBarostat()) |
1587 |
|
myNPTf->setTauBarostat(globals->getTauBarostat()); |
1462 |
– |
else{ |
1463 |
– |
sprintf(painCave.errMsg, |
1464 |
– |
"SimSetup error: If you use an NPT\n" |
1465 |
– |
" ensemble, you must set tauBarostat.\n"); |
1466 |
– |
painCave.isFatal = 1; |
1467 |
– |
simError(); |
1468 |
– |
} |
1469 |
– |
break; |
1588 |
|
|
1471 |
– |
case NPTim_ENS: |
1472 |
– |
if (globals->haveZconstraints()){ |
1473 |
– |
setupZConstraint(info[k]); |
1474 |
– |
myNPTim = new ZConstraint<NPTim<RealIntegrator> >(&(info[k]), the_ff); |
1475 |
– |
} |
1476 |
– |
else |
1477 |
– |
myNPTim = new NPTim<RealIntegrator>(&(info[k]), the_ff); |
1478 |
– |
|
1479 |
– |
myNPTim->setTargetTemp(globals->getTargetTemp()); |
1480 |
– |
|
1481 |
– |
if (globals->haveTargetPressure()) |
1482 |
– |
myNPTim->setTargetPressure(globals->getTargetPressure()); |
1589 |
|
else{ |
1590 |
|
sprintf(painCave.errMsg, |
1485 |
– |
"SimSetup error: If you use a constant pressure\n" |
1486 |
– |
" ensemble, you must set targetPressure in the BASS file.\n"); |
1487 |
– |
painCave.isFatal = 1; |
1488 |
– |
simError(); |
1489 |
– |
} |
1490 |
– |
|
1491 |
– |
if (globals->haveTauThermostat()) |
1492 |
– |
myNPTim->setTauThermostat(globals->getTauThermostat()); |
1493 |
– |
else{ |
1494 |
– |
sprintf(painCave.errMsg, |
1591 |
|
"SimSetup error: If you use an NPT\n" |
1592 |
< |
" ensemble, you must set tauThermostat.\n"); |
1592 |
> |
"\tensemble, you must set tauBarostat.\n"); |
1593 |
|
painCave.isFatal = 1; |
1594 |
|
simError(); |
1595 |
|
} |
1596 |
|
|
1597 |
< |
if (globals->haveTauBarostat()) |
1502 |
< |
myNPTim->setTauBarostat(globals->getTauBarostat()); |
1503 |
< |
else{ |
1504 |
< |
sprintf(painCave.errMsg, |
1505 |
< |
"SimSetup error: If you use an NPT\n" |
1506 |
< |
" ensemble, you must set tauBarostat.\n"); |
1507 |
< |
painCave.isFatal = 1; |
1508 |
< |
simError(); |
1509 |
< |
} |
1597 |
> |
info->the_integrator = myNPTf; |
1598 |
|
break; |
1599 |
|
|
1600 |
< |
case NPTfm_ENS: |
1600 |
> |
case NPTxyz_ENS: |
1601 |
|
if (globals->haveZconstraints()){ |
1602 |
|
setupZConstraint(info[k]); |
1603 |
< |
myNPTfm = new ZConstraint<NPTfm<RealIntegrator> >(&(info[k]), the_ff); |
1603 |
> |
myNPTxyz = new ZConstraint<NPTxyz<NPT <RealIntegrator> > >(&(info[k]), the_ff); |
1604 |
|
} |
1605 |
|
else |
1606 |
< |
myNPTfm = new NPTfm<RealIntegrator>(&(info[k]), the_ff); |
1606 |
> |
myNPTxyz = new NPTxyz<NPT <RealIntegrator> >(&(info[k]), the_ff); |
1607 |
|
|
1608 |
< |
myNPTfm->setTargetTemp(globals->getTargetTemp()); |
1608 |
> |
myNPTxyz->setTargetTemp(globals->getTargetTemp()); |
1609 |
|
|
1610 |
|
if (globals->haveTargetPressure()) |
1611 |
< |
myNPTfm->setTargetPressure(globals->getTargetPressure()); |
1611 |
> |
myNPTxyz->setTargetPressure(globals->getTargetPressure()); |
1612 |
|
else{ |
1613 |
|
sprintf(painCave.errMsg, |
1614 |
|
"SimSetup error: If you use a constant pressure\n" |
1615 |
< |
" ensemble, you must set targetPressure in the BASS file.\n"); |
1615 |
> |
"\tensemble, you must set targetPressure in the BASS file.\n"); |
1616 |
|
painCave.isFatal = 1; |
1617 |
|
simError(); |
1618 |
< |
} |
1618 |
> |
} |
1619 |
|
|
1620 |
|
if (globals->haveTauThermostat()) |
1621 |
< |
myNPTfm->setTauThermostat(globals->getTauThermostat()); |
1621 |
> |
myNPTxyz->setTauThermostat(globals->getTauThermostat()); |
1622 |
|
else{ |
1623 |
|
sprintf(painCave.errMsg, |
1624 |
|
"SimSetup error: If you use an NPT\n" |
1625 |
< |
" ensemble, you must set tauThermostat.\n"); |
1625 |
> |
"\tensemble, you must set tauThermostat.\n"); |
1626 |
|
painCave.isFatal = 1; |
1627 |
|
simError(); |
1628 |
|
} |
1629 |
|
|
1630 |
|
if (globals->haveTauBarostat()) |
1631 |
< |
myNPTfm->setTauBarostat(globals->getTauBarostat()); |
1631 |
> |
myNPTxyz->setTauBarostat(globals->getTauBarostat()); |
1632 |
|
else{ |
1633 |
|
sprintf(painCave.errMsg, |
1634 |
|
"SimSetup error: If you use an NPT\n" |
1635 |
< |
" ensemble, you must set tauBarostat.\n"); |
1635 |
> |
"\tensemble, you must set tauBarostat.\n"); |
1636 |
|
painCave.isFatal = 1; |
1637 |
|
simError(); |
1638 |
|
} |
1639 |
+ |
|
1640 |
+ |
info->the_integrator = myNPTxyz; |
1641 |
|
break; |
1642 |
|
|
1643 |
|
default: |
1685 |
|
} |
1686 |
|
else{ |
1687 |
|
sprintf(painCave.errMsg, |
1688 |
< |
"ZConstraint error: If you use an ZConstraint\n" |
1689 |
< |
" , you must set sample time.\n"); |
1688 |
> |
"ZConstraint error: If you use a ZConstraint,\n" |
1689 |
> |
"\tyou must set zconsTime.\n"); |
1690 |
|
painCave.isFatal = 1; |
1691 |
|
simError(); |
1692 |
|
} |
1701 |
|
else{ |
1702 |
|
double defaultZConsTol = 0.01; |
1703 |
|
sprintf(painCave.errMsg, |
1704 |
< |
"ZConstraint Waring: Tolerance for z-constraint methodl is not specified\n" |
1705 |
< |
" , default value %f is used.\n", |
1704 |
> |
"ZConstraint Warning: Tolerance for z-constraint method is not specified.\n" |
1705 |
> |
"\tOOPSE will use a default value of %f.\n" |
1706 |
> |
"\tTo set the tolerance, use the zconsTol variable.\n", |
1707 |
|
defaultZConsTol); |
1708 |
|
painCave.isFatal = 0; |
1709 |
|
simError(); |
1721 |
|
} |
1722 |
|
else{ |
1723 |
|
sprintf(painCave.errMsg, |
1724 |
< |
"ZConstraint Warning: User does not set force Subtraction policy, " |
1725 |
< |
"PolicyByMass is used\n"); |
1724 |
> |
"ZConstraint Warning: No force subtraction policy was set.\n" |
1725 |
> |
"\tOOPSE will use PolicyByMass.\n" |
1726 |
> |
"\tTo set the policy, use the zconsForcePolicy variable.\n"); |
1727 |
|
painCave.isFatal = 0; |
1728 |
|
simError(); |
1729 |
|
zconsForcePolicy->setData("BYMASS"); |
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 |
|
|
1795 |
|
//check the uniqueness of index |
1796 |
|
if(!zconsParaData->isIndexUnique()){ |
1797 |
|
sprintf(painCave.errMsg, |
1798 |
< |
"ZConstraint Error: molIndex is not unique\n"); |
1798 |
> |
"ZConstraint Error: molIndex is not unique!\n"); |
1799 |
|
painCave.isFatal = 1; |
1800 |
|
simError(); |
1801 |
|
} |
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 |
+ |
} |