ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/SimSetup.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 1035 by tim, Fri Feb 6 21:37:59 2004 UTC vs.
Revision 1116 by tim, Thu Apr 15 22:15:21 2004 UTC

# Line 9 | Line 9
9   #include "parse_me.h"
10   #include "Integrator.hpp"
11   #include "simError.h"
12 < #include "ConjugateMinimizer.hpp"
12 > #include "RigidBody.hpp"
13 > //#include "ConjugateMinimizer.hpp"
14 > #include "OOPSEMinimizer.hpp"
15  
16   #ifdef IS_MPI
17   #include "mpiBASS.h"
# Line 28 | Line 30
30   #define FF_DUFF  0
31   #define FF_LJ    1
32   #define FF_EAM   2
33 < #define FF_H2O 3
33 > #define FF_H2O   3
34  
35   using namespace std;
36  
# Line 145 | Line 147 | void SimSetup::createSim(void){
147    // make the output filenames
148  
149    makeOutNames();
148 –
149 –  if (globals->haveMinimizer())
150 –    // make minimizer
151 –    makeMinimizer();
152 –  else
153 –    // make the integrator
154 –    makeIntegrator();
150    
151   #ifdef IS_MPI
152    mpiSim->mpiRefresh();
# Line 160 | Line 155 | void SimSetup::createSim(void){
155    // initialize the Fortran
156  
157    initFortran();
158 +
159 +  if (globals->haveMinimizer())
160 +    // make minimizer
161 +    makeMinimizer();
162 +  else
163 +    // make the integrator
164 +    makeIntegrator();
165 +
166   }
167  
168  
169   void SimSetup::makeMolecules(void){
170 <  int k;
171 <  int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset;
170 >  int i, j, k;
171 >  int exI, exJ, exK, exL, slI, slJ;
172 >  int tempI, tempJ, tempK, tempL;
173 >  int molI;
174 >  int stampID, atomOffset, rbOffset;
175    molInit molInfo;
176    DirectionalAtom* dAtom;
177 +  RigidBody* myRB;
178 +  StuntDouble* mySD;
179    LinkedAssign* extras;
180    LinkedAssign* current_extra;
181    AtomStamp* currentAtom;
182    BondStamp* currentBond;
183    BendStamp* currentBend;
184    TorsionStamp* currentTorsion;
185 +  RigidBodyStamp* currentRigidBody;
186  
187    bond_pair* theBonds;
188    bend_set* theBends;
189    torsion_set* theTorsions;
190  
191 +  set<int> skipList;
192 +
193 +  double phi, theta, psi;
194 +  char* molName;
195 +  char rbName[100];
196 +
197    //init the forceField paramters
198  
199    the_ff->readParams();
200  
186 –
201    // init the atoms
202  
203 <  double phi, theta, psi;
190 <  double sux, suy, suz;
191 <  double Axx, Axy, Axz, Ayx, Ayy, Ayz, Azx, Azy, Azz;
192 <  double ux, uy, uz, u, uSqr;
203 >  int nMembers, nNew, rb1, rb2;
204  
205    for (k = 0; k < nInfo; k++){
206      the_ff->setSimInfo(&(info[k]));
207  
208      atomOffset = 0;
209 <    excludeOffset = 0;
209 >
210      for (i = 0; i < info[k].n_mol; i++){
211        stampID = info[k].molecules[i].getStampID();
212 +      molName = comp_stamps[stampID]->getID();
213  
214        molInfo.nAtoms = comp_stamps[stampID]->getNAtoms();
215        molInfo.nBonds = comp_stamps[stampID]->getNBonds();
216        molInfo.nBends = comp_stamps[stampID]->getNBends();
217        molInfo.nTorsions = comp_stamps[stampID]->getNTorsions();
218 <      molInfo.nExcludes = molInfo.nBonds + molInfo.nBends + molInfo.nTorsions;
219 <
218 >      molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies();
219 >      
220        molInfo.myAtoms = &(info[k].atoms[atomOffset]);
209 –      molInfo.myExcludes = &(info[k].excludes[excludeOffset]);
210 –      molInfo.myBonds = new Bond * [molInfo.nBonds];
211 –      molInfo.myBends = new Bend * [molInfo.nBends];
212 –      molInfo.myTorsions = new Torsion * [molInfo.nTorsions];
221  
222 +      if (molInfo.nBonds > 0)
223 +        molInfo.myBonds = new (Bond *) [molInfo.nBonds];
224 +      else
225 +        molInfo.myBonds = NULL;
226 +
227 +      if (molInfo.nBends > 0)
228 +        molInfo.myBends = new (Bend *) [molInfo.nBends];
229 +      else
230 +        molInfo.myBends = NULL;
231 +
232 +      if (molInfo.nTorsions > 0)
233 +        molInfo.myTorsions = new (Torsion *) [molInfo.nTorsions];
234 +      else
235 +        molInfo.myTorsions = NULL;
236 +
237        theBonds = new bond_pair[molInfo.nBonds];
238        theBends = new bend_set[molInfo.nBends];
239        theTorsions = new torsion_set[molInfo.nTorsions];
240 <
240 >      
241        // make the Atoms
242  
243        for (j = 0; j < molInfo.nAtoms; j++){
244          currentAtom = comp_stamps[stampID]->getAtom(j);
245 +
246          if (currentAtom->haveOrientation()){
247            dAtom = new DirectionalAtom((j + atomOffset),
248                                        info[k].getConfiguration());
# Line 232 | Line 256 | void SimSetup::makeMolecules(void){
256            phi = currentAtom->getEulerPhi() * M_PI / 180.0;
257            theta = currentAtom->getEulerTheta() * M_PI / 180.0;
258            psi = currentAtom->getEulerPsi()* M_PI / 180.0;
259 +
260 +          dAtom->setUnitFrameFromEuler(phi, theta, psi);
261              
262 <          Axx = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
263 <          Axy = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
238 <          Axz = sin(theta) * sin(psi);
239 <          
240 <          Ayx = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
241 <          Ayy = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
242 <          Ayz = sin(theta) * cos(psi);
243 <          
244 <          Azx = sin(phi) * sin(theta);
245 <          Azy = -cos(phi) * sin(theta);
246 <          Azz = cos(theta);
262 >        }
263 >        else{
264  
265 <          sux = 0.0;
249 <          suy = 0.0;
250 <          suz = 1.0;
265 >          molInfo.myAtoms[j] = new Atom((j + atomOffset), info[k].getConfiguration());
266  
252 –          ux = (Axx * sux) + (Ayx * suy) + (Azx * suz);
253 –          uy = (Axy * sux) + (Ayy * suy) + (Azy * suz);
254 –          uz = (Axz * sux) + (Ayz * suy) + (Azz * suz);
255 –
256 –          uSqr = (ux * ux) + (uy * uy) + (uz * uz);
257 –
258 –          u = sqrt(uSqr);
259 –          ux = ux / u;
260 –          uy = uy / u;
261 –          uz = uz / u;
262 –
263 –          dAtom->setSUx(ux);
264 –          dAtom->setSUy(uy);
265 –          dAtom->setSUz(uz);
267          }
267 –        else{
268 –          molInfo.myAtoms[j] = new GeneralAtom((j + atomOffset),
269 –                                               info[k].getConfiguration());
270 –        }
271 –        molInfo.myAtoms[j]->setType(currentAtom->getType());
268  
269 +        molInfo.myAtoms[j]->setType(currentAtom->getType());
270   #ifdef IS_MPI
271  
272 <        molInfo.myAtoms[j]->setGlobalIndex(globalIndex[j + atomOffset]);
272 >        molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]);
273  
274   #endif // is_mpi
275        }
# Line 283 | Line 280 | void SimSetup::makeMolecules(void){
280          theBonds[j].a = currentBond->getA() + atomOffset;
281          theBonds[j].b = currentBond->getB() + atomOffset;
282  
283 <        exI = theBonds[j].a;
284 <        exJ = theBonds[j].b;
283 >        tempI = theBonds[j].a;
284 >        tempJ = theBonds[j].b;
285  
289 –        // exclude_I must always be the smaller of the pair
290 –        if (exI > exJ){
291 –          tempEx = exI;
292 –          exI = exJ;
293 –          exJ = tempEx;
294 –        }
286   #ifdef IS_MPI
287 <        tempEx = exI;
288 <        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
289 <        tempEx = exJ;
290 <        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
291 <
292 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
302 < #else  // isn't MPI
287 >        exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
288 >        exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
289 > #else
290 >        exI = tempI + 1;
291 >        exJ = tempJ + 1;
292 > #endif
293  
294 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
305 < #endif  //is_mpi
294 >        info[k].excludes->addPair(exI, exJ);
295        }
307 –      excludeOffset += molInfo.nBonds;
296  
297        //make the bends
298        for (j = 0; j < molInfo.nBends; j++){
# Line 354 | Line 342 | void SimSetup::makeMolecules(void){
342            }
343          }
344  
345 <        if (!theBends[j].isGhost){
346 <          exI = theBends[j].a;
347 <          exJ = theBends[j].c;
348 <        }
349 <        else{
362 <          exI = theBends[j].a;
363 <          exJ = theBends[j].b;
364 <        }
365 <
366 <        // exclude_I must always be the smaller of the pair
367 <        if (exI > exJ){
368 <          tempEx = exI;
369 <          exI = exJ;
370 <          exJ = tempEx;
371 <        }
345 >        if (theBends[j].isGhost) {
346 >          
347 >          tempI = theBends[j].a;
348 >          tempJ = theBends[j].b;
349 >          
350   #ifdef IS_MPI
351 <        tempEx = exI;
352 <        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
353 <        tempEx = exJ;
354 <        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
351 >          exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
352 >          exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
353 > #else
354 >          exI = tempI + 1;
355 >          exJ = tempJ + 1;
356 > #endif          
357 >          info[k].excludes->addPair(exI, exJ);
358  
359 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
360 < #else  // isn't MPI
361 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
362 < #endif  //is_mpi
359 >        } else {
360 >
361 >          tempI = theBends[j].a;
362 >          tempJ = theBends[j].b;
363 >          tempK = theBends[j].c;
364 >          
365 > #ifdef IS_MPI
366 >          exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
367 >          exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
368 >          exK = info[k].atoms[tempK]->getGlobalIndex() + 1;
369 > #else
370 >          exI = tempI + 1;
371 >          exJ = tempJ + 1;
372 >          exK = tempK + 1;
373 > #endif
374 >          
375 >          info[k].excludes->addPair(exI, exK);
376 >          info[k].excludes->addPair(exI, exJ);
377 >          info[k].excludes->addPair(exJ, exK);
378 >        }
379        }
383 –      excludeOffset += molInfo.nBends;
380  
381        for (j = 0; j < molInfo.nTorsions; j++){
382          currentTorsion = comp_stamps[stampID]->getTorsion(j);
# Line 389 | Line 385 | void SimSetup::makeMolecules(void){
385          theTorsions[j].c = currentTorsion->getC() + atomOffset;
386          theTorsions[j].d = currentTorsion->getD() + atomOffset;
387  
388 <        exI = theTorsions[j].a;
389 <        exJ = theTorsions[j].d;
388 >        tempI = theTorsions[j].a;      
389 >        tempJ = theTorsions[j].b;
390 >        tempK = theTorsions[j].c;
391 >        tempL = theTorsions[j].d;
392  
395 –        // exclude_I must always be the smaller of the pair
396 –        if (exI > exJ){
397 –          tempEx = exI;
398 –          exI = exJ;
399 –          exJ = tempEx;
400 –        }
393   #ifdef IS_MPI
394 <        tempEx = exI;
395 <        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
396 <        tempEx = exJ;
397 <        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
394 >        exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
395 >        exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
396 >        exK = info[k].atoms[tempK]->getGlobalIndex() + 1;
397 >        exL = info[k].atoms[tempL]->getGlobalIndex() + 1;
398 > #else
399 >        exI = tempI + 1;
400 >        exJ = tempJ + 1;
401 >        exK = tempK + 1;
402 >        exL = tempL + 1;
403 > #endif
404  
405 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
406 < #else  // isn't MPI
407 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
408 < #endif  //is_mpi
405 >        info[k].excludes->addPair(exI, exJ);
406 >        info[k].excludes->addPair(exI, exK);
407 >        info[k].excludes->addPair(exI, exL);        
408 >        info[k].excludes->addPair(exJ, exK);
409 >        info[k].excludes->addPair(exJ, exL);
410 >        info[k].excludes->addPair(exK, exL);
411        }
412 –      excludeOffset += molInfo.nTorsions;
412  
413 +      
414 +      molInfo.myRigidBodies.clear();
415 +      
416 +      for (j = 0; j < molInfo.nRigidBodies; j++){
417  
418 <      // send the arrays off to the forceField for init.
418 >        currentRigidBody = comp_stamps[stampID]->getRigidBody(j);
419 >        nMembers = currentRigidBody->getNMembers();
420 >
421 >        // Create the Rigid Body:
422 >
423 >        myRB = new RigidBody();
424 >
425 >        sprintf(rbName,"%s_RB_%d", molName, j);
426 >        myRB->setType(rbName);
427 >        
428 >        for (rb1 = 0; rb1 < nMembers; rb1++) {
429 >
430 >          // molI is atom numbering inside this molecule
431 >          molI = currentRigidBody->getMember(rb1);    
432  
433 +          // tempI is atom numbering on local processor
434 +          tempI = molI + atomOffset;
435 +
436 +          // currentAtom is the AtomStamp (which we need for
437 +          // rigid body reference positions)
438 +          currentAtom = comp_stamps[stampID]->getAtom(molI);
439 +
440 +          // When we add to the rigid body, add the atom itself and
441 +          // the stamp info:
442 +
443 +          myRB->addAtom(info[k].atoms[tempI], currentAtom);
444 +          
445 +          // Add this atom to the Skip List for the integrators
446 + #ifdef IS_MPI
447 +          slI = info[k].atoms[tempI]->getGlobalIndex();
448 + #else
449 +          slI = tempI;
450 + #endif
451 +          skipList.insert(slI);
452 +          
453 +        }
454 +        
455 +        for(rb1 = 0; rb1 < nMembers - 1; rb1++) {
456 +          for(rb2 = rb1+1; rb2 < nMembers; rb2++) {
457 +            
458 +            tempI = currentRigidBody->getMember(rb1);
459 +            tempJ = currentRigidBody->getMember(rb2);
460 +            
461 +            // Some explanation is required here.
462 +            // Fortran indexing starts at 1, while c indexing starts at 0
463 +            // Also, in parallel computations, the GlobalIndex is
464 +            // used for the exclude list:
465 +            
466 + #ifdef IS_MPI
467 +            exI = molInfo.myAtoms[tempI]->getGlobalIndex() + 1;
468 +            exJ = molInfo.myAtoms[tempJ]->getGlobalIndex() + 1;
469 + #else
470 +            exI = molInfo.myAtoms[tempI]->getIndex() + 1;
471 +            exJ = molInfo.myAtoms[tempJ]->getIndex() + 1;
472 + #endif
473 +            
474 +            info[k].excludes->addPair(exI, exJ);
475 +            
476 +          }
477 +        }
478 +
479 +        molInfo.myRigidBodies.push_back(myRB);
480 +        info[k].rigidBodies.push_back(myRB);
481 +      }
482 +      
483 +
484 +      // After this is all set up, scan through the atoms to
485 +      // see if they can be added to the integrableObjects:
486 +
487 +      molInfo.myIntegrableObjects.clear();
488 +      
489 +
490 +      for (j = 0; j < molInfo.nAtoms; j++){
491 +
492 + #ifdef IS_MPI
493 +        slJ = molInfo.myAtoms[j]->getGlobalIndex();
494 + #else
495 +        slJ = j+atomOffset;
496 + #endif
497 +
498 +        // if they aren't on the skip list, then they can be integrated
499 +
500 +        if (skipList.find(slJ) == skipList.end()) {
501 +          mySD = (StuntDouble *) molInfo.myAtoms[j];
502 +          info[k].integrableObjects.push_back(mySD);
503 +          molInfo.myIntegrableObjects.push_back(mySD);
504 +        }
505 +      }
506 +
507 +      // all rigid bodies are integrated:
508 +
509 +      for (j = 0; j < molInfo.nRigidBodies; j++) {
510 +        mySD = (StuntDouble *) molInfo.myRigidBodies[j];
511 +        info[k].integrableObjects.push_back(mySD);      
512 +        molInfo.myIntegrableObjects.push_back(mySD);
513 +      }
514 +    
515 +      
516 +      // send the arrays off to the forceField for init.
517 +      
518        the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms);
519        the_ff->initializeBonds(molInfo.nBonds, molInfo.myBonds, theBonds);
520        the_ff->initializeBends(molInfo.nBends, molInfo.myBends, theBends);
521        the_ff->initializeTorsions(molInfo.nTorsions, molInfo.myTorsions,
522                                   theTorsions);
523  
423 –
524        info[k].molecules[i].initialize(molInfo);
525  
526  
# Line 428 | Line 528 | void SimSetup::makeMolecules(void){
528        delete[] theBonds;
529        delete[] theBends;
530        delete[] theTorsions;
531 <    }
531 >    }    
532    }
533  
534   #ifdef IS_MPI
# Line 438 | Line 538 | void SimSetup::makeMolecules(void){
538  
539    // clean up the forcefield
540  
541 <  the_ff->calcRcut();
541 >  if (!globals->haveLJrcut()){
542 >
543 >    the_ff->calcRcut();
544 >
545 >  } else {
546 >    
547 >    the_ff->setRcut( globals->getLJrcut() );
548 >  }
549 >
550    the_ff->cleanMe();
551   }
552  
# Line 856 | Line 964 | void SimSetup::finalInfoCheck(void){
964   void SimSetup::finalInfoCheck(void){
965    int index;
966    int usesDipoles;
967 +  int usesCharges;
968    int i;
969  
970    for (i = 0; i < nInfo; i++){
# Line 867 | Line 976 | void SimSetup::finalInfoCheck(void){
976        usesDipoles = (info[i].atoms[index])->hasDipole();
977        index++;
978      }
979 <
979 >    index = 0;
980 >    usesCharges = 0;
981 >    while ((index < info[i].n_atoms) && !usesCharges){
982 >      usesCharges= (info[i].atoms[index])->hasCharge();
983 >      index++;
984 >    }
985   #ifdef IS_MPI
986      int myUse = usesDipoles;
987      MPI_Allreduce(&myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
# Line 918 | Line 1032 | void SimSetup::finalInfoCheck(void){
1032        info[i].dielectric = globals->getDielectric();
1033      }
1034      else{
1035 <      if (usesDipoles){
1035 >      if (usesDipoles || usesCharges){
1036          if (!globals->haveECR()){
1037            sprintf(painCave.errMsg,
1038                    "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
# Line 1220 | Line 1334 | void SimSetup::calcSysValues(void){
1334    tot_bonds = 0;
1335    tot_bends = 0;
1336    tot_torsions = 0;
1337 +  tot_rigid = 0;
1338    for (i = 0; i < n_components; i++){
1339      tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
1340      tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
1341      tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
1342      tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
1343 +    tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies();
1344    }
1345 <
1345 >  
1346    tot_SRI = tot_bonds + tot_bends + tot_torsions;
1347    molMembershipArray = new int[tot_atoms];
1348  
# Line 1248 | Line 1364 | void SimSetup::mpiMolDivide(void){
1364    int i, j, k;
1365    int localMol, allMol;
1366    int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1367 +  int local_rigid;
1368 +  vector<int> globalAtomIndex;
1369 +  vector<int> globalMolIndex;
1370  
1371    mpiSim = new mpiSimulation(info);
1372  
1373 <  globalIndex = mpiSim->divideLabor();
1373 >  mpiSim->divideLabor();
1374 >  globalAtomIndex = mpiSim->getGlobalAtomIndex();
1375 >  globalMolIndex = mpiSim->getGlobalMolIndex();
1376  
1377    // set up the local variables
1378  
# Line 1264 | Line 1385 | void SimSetup::mpiMolDivide(void){
1385    local_bonds = 0;
1386    local_bends = 0;
1387    local_torsions = 0;
1388 <  globalAtomIndex = 0;
1388 >  local_rigid = 0;
1389 >  globalAtomCounter = 0;
1390  
1269 –
1391    for (i = 0; i < n_components; i++){
1392      for (j = 0; j < components_nmol[i]; j++){
1393        if (mol2proc[allMol] == worldRank){
# Line 1274 | Line 1395 | void SimSetup::mpiMolDivide(void){
1395          local_bonds += comp_stamps[i]->getNBonds();
1396          local_bends += comp_stamps[i]->getNBends();
1397          local_torsions += comp_stamps[i]->getNTorsions();
1398 +        local_rigid += comp_stamps[i]->getNRigidBodies();
1399          localMol++;
1400        }      
1401        for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1402 <        info[0].molMembershipArray[globalAtomIndex] = allMol;
1403 <        globalAtomIndex++;
1402 >        info[0].molMembershipArray[globalAtomCounter] = allMol;
1403 >        globalAtomCounter++;
1404        }
1405  
1406        allMol++;
# Line 1287 | Line 1409 | void SimSetup::mpiMolDivide(void){
1409    local_SRI = local_bonds + local_bends + local_torsions;
1410  
1411    info[0].n_atoms = mpiSim->getMyNlocal();  
1412 +  
1413  
1414    if (local_atoms != info[0].n_atoms){
1415      sprintf(painCave.errMsg,
# Line 1319 | Line 1442 | void SimSetup::makeSysArrays(void){
1442  
1443    Atom** the_atoms;
1444    Molecule* the_molecules;
1322 –  Exclude** the_excludes;
1445  
1324 –
1446    for (l = 0; l < nInfo; l++){
1447      // create the atom and short range interaction arrays
1448  
# Line 1347 | Line 1468 | void SimSetup::makeSysArrays(void){
1468   #else // is_mpi
1469  
1470      molIndex = 0;
1471 <    globalAtomIndex = 0;
1471 >    globalAtomCounter = 0;
1472      for (i = 0; i < n_components; i++){
1473        for (j = 0; j < components_nmol[i]; j++){
1474          the_molecules[molIndex].setStampID(i);
1475          the_molecules[molIndex].setMyIndex(molIndex);
1476          the_molecules[molIndex].setGlobalIndex(molIndex);
1477          for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1478 <          info[l].molMembershipArray[globalAtomIndex] = molIndex;
1479 <          globalAtomIndex++;
1478 >          info[l].molMembershipArray[globalAtomCounter] = molIndex;
1479 >          globalAtomCounter++;
1480          }
1481          molIndex++;
1482        }
# Line 1363 | Line 1484 | void SimSetup::makeSysArrays(void){
1484  
1485  
1486   #endif // is_mpi
1366 –
1367 –
1368 –    if (info[l].n_SRI){
1369 –      Exclude::createArray(info[l].n_SRI);
1370 –      the_excludes = new Exclude * [info[l].n_SRI];
1371 –      for (int ex = 0; ex < info[l].n_SRI; ex++){
1372 –        the_excludes[ex] = new Exclude(ex);
1373 –      }
1374 –      info[l].globalExcludes = new int;
1375 –      info[l].n_exclude = info[l].n_SRI;
1376 –    }
1377 –    else{
1378 –      Exclude::createArray(1);
1379 –      the_excludes = new Exclude * ;
1380 –      the_excludes[0] = new Exclude(0);
1381 –      the_excludes[0]->setPair(0, 0);
1382 –      info[l].globalExcludes = new int;
1383 –      info[l].globalExcludes[0] = 0;
1384 –      info[l].n_exclude = 0;
1385 –    }
1487  
1488 +    info[l].globalExcludes = new int;
1489 +    info[l].globalExcludes[0] = 0;
1490 +    
1491      // set the arrays into the SimInfo object
1492  
1493      info[l].atoms = the_atoms;
1494      info[l].molecules = the_molecules;
1495      info[l].nGlobalExcludes = 0;
1392 –    info[l].excludes = the_excludes;
1496  
1497      the_ff->setSimInfo(info);
1498    }
# Line 1662 | Line 1765 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1765    }
1766  
1767    theInfo.addProperty(zconsForcePolicy);
1768 +
1769 +  //set zcons gap
1770 +  DoubleData* zconsGap = new DoubleData();
1771 +  zconsGap->setID(ZCONSGAP_ID);
1772 +
1773 +  if (globals->haveZConsGap()){
1774 +    zconsGap->setData(globals->getZconsGap());
1775 +    theInfo.addProperty(zconsGap);  
1776 +  }
1777 +
1778 +  //set zcons fixtime
1779 +  DoubleData* zconsFixtime = new DoubleData();
1780 +  zconsFixtime->setID(ZCONSFIXTIME_ID);
1781 +
1782 +  if (globals->haveZConsFixTime()){
1783 +    zconsFixtime->setData(globals->getZconsFixtime());
1784 +    theInfo.addProperty(zconsFixtime);  
1785 +  }
1786 +
1787 +  //set zconsUsingSMD
1788 +  IntData* zconsUsingSMD = new IntData();
1789 +  zconsUsingSMD->setID(ZCONSUSINGSMD_ID);
1790 +
1791 +  if (globals->haveZConsUsingSMD()){
1792 +    zconsUsingSMD->setData(globals->getZconsUsingSMD());
1793 +    theInfo.addProperty(zconsUsingSMD);  
1794 +  }
1795  
1796    //Determine the name of ouput file and add it into SimInfo's property list
1797    //Be careful, do not use inFileName, since it is a pointer which
# Line 1692 | Line 1822 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1822      tempParaItem.zPos = zconStamp[i]->getZpos();
1823      tempParaItem.zconsIndex = zconStamp[i]->getMolIndex();
1824      tempParaItem.kRatio = zconStamp[i]->getKratio();
1825 <
1825 >    tempParaItem.havingCantVel = zconStamp[i]->haveCantVel();
1826 >    tempParaItem.cantVel = zconStamp[i]->getCantVel();    
1827      zconsParaData->addItem(tempParaItem);
1828    }
1829  
# Line 1713 | Line 1844 | void SimSetup::makeMinimizer(){
1844  
1845   void SimSetup::makeMinimizer(){
1846  
1847 <  OOPSEMinimizerBase* myOOPSEMinimizerBase;
1717 <  ObjFunctor1 * objFunc;
1718 <  OutputFunctor* outputFunc;
1719 <  ConcreteNLModel1* nlp;
1847 >  OOPSEMinimizer* myOOPSEMinimizer;
1848    MinimizerParameterSet* param;
1849 <  ConjugateMinimizerBase* minimizer;
1722 <  int dim;
1849 >  char minimizerName[100];
1850    
1851    for (int i = 0; i < nInfo; i++){
1852 <    //creat
1726 <    myOOPSEMinimizerBase = new OOPSEMinimizerBase(&(info[i]), the_ff);
1727 <
1728 <     info[i].the_integrator = myOOPSEMinimizerBase;
1729 <    //creat the object functor;
1730 <    objFunc = (ObjFunctor1*) new ClassMemObjFunctor1<OOPSEMinimizerBase>
1731 <                                              (myOOPSEMinimizerBase, &OOPSEMinimizerBase::calcGradient);
1732 <
1733 <    //creat output functor;
1734 <    outputFunc =  new ClassMemOutputFunctor<OOPSEMinimizerBase>
1735 <                               (myOOPSEMinimizerBase, &OOPSEMinimizerBase::output);
1736 <
1737 <    //creat nonlinear model
1738 <    dim = myOOPSEMinimizerBase->getDim();    
1739 <    nlp = new ConcreteNLModel1(dim, objFunc);
1740 <
1741 <    nlp->setX(myOOPSEMinimizerBase->getCoor());
1742 <
1852 >    
1853      //prepare parameter set for minimizer
1854      param = new MinimizerParameterSet();
1855      param->setDefaultParameter();
# Line 1764 | Line 1874 | void SimSetup::makeMinimizer(){
1874        param->setWriteFrq(globals->getMinWriteFrq());
1875      }
1876      
1877 <    if (globals->haveMinResetFrq()){
1878 <      param->setResetFrq(globals->getMinResetFrq());
1877 >    if (globals->haveMinStepSize()){
1878 >      param->setStepSize(globals->getMinStepSize());
1879      }
1880  
1881      if (globals->haveMinLSMaxIter()){
# Line 1775 | Line 1885 | void SimSetup::makeMinimizer(){
1885      if (globals->haveMinLSTol()){
1886        param->setLineSearchTol(globals->getMinLSTol());
1887      }    
1778 –    
1779 –     //creat the minimizer
1780 –     minimizer = new PRCGMinimizer(nlp, param);
1781 –     minimizer->setLineSearchStrategy(nlp, GoldenSection);
1782 –     minimizer->setOutputFunctor(outputFunc);
1888  
1889 +    strcpy(minimizerName, globals->getMinimizer());
1890 +
1891 +    if (!strcasecmp(minimizerName, "CG")){
1892 +      myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param);
1893 +    }
1894 +    else if (!strcasecmp(minimizerName, "SD")){
1895 +    //myOOPSEMinimizer = MinimizerFactory.creatMinimizer("", &(info[i]), the_ff, param);
1896 +      myOOPSEMinimizer = new SDMinimizer(&(info[i]), the_ff, param);
1897 +    }
1898 +    else{
1899 +          sprintf(painCave.errMsg,
1900 +                  "SimSetup error: Unrecognized Minimizer, use Conjugate Gradient \n");
1901 +          painCave.isFatal = 0;
1902 +          simError();
1903 +
1904 +      myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param);          
1905 +    }
1906 +     info[i].the_integrator = myOOPSEMinimizer;
1907 +
1908       //store the minimizer into simInfo
1909 <     info[i].the_minimizer = minimizer;
1909 >     info[i].the_minimizer = myOOPSEMinimizer;
1910       info[i].has_minimizer = true;
1911    }
1912  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines