[ViewVC] Diff of: group/trunk/OOPSE/libmdtools/SimSetup.cpp

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 780 by mmeineke, Mon Sep 22 21:23:25 2003 UTC vs.
Revision 1097 by gezelter, Mon Apr 12 20:32:20 2004 UTC

#	Line 1 \| Line 1
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"
#	Line 22 \| Line 24
24		#define NVT_ENS 1
25		#define NPTi_ENS 2
26		#define NPTf_ENS 3
27	+	#define NPTxyz_ENS 4
28
26	–	#define FF_DUFF 0
27	–	#define FF_LJ 1
28	–	#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
#	Line 52 \| Line 83 \| *void SimSetup::setSimInfo(SimInfo the_info, int theNi**
83		info = the_info;
84		nInfo = theNinfo;
85		isInfoArray = 1;
86	+	initSuspend = true;
87		}
88
89
#	Line 90 \| Line 122 \| void SimSetup::createSim(void){
122		#endif // is_mpi
123
124		void SimSetup::createSim(void){
93	–	int i, j, k, globalAtomIndex;
125
126		// gather all of the information from the Bass file
127
#	Line 106 \| Line 137 \| void SimSetup::createSim(void){
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
#	Line 129 \| Line 166 \| void SimSetup::makeMolecules(void){
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
152	–
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
#	Line 166 \| Line 211 \| void SimSetup::makeMolecules(void){
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]);
172	–	molInfo.myExcludes = &(info[k].excludes[excludeOffset]);
173	–	molInfo.myBonds = new Bond * [molInfo.nBonds];
174	–	molInfo.myBends = new Bend * [molInfo.nBends];
175	–	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();
194	<
195	<	uSqr = (ux * ux) + (uy * uy) + (uz * uz);
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	<	u = sqrt(uSqr);
253	<	ux = ux / u;
254	<	uy = uy / u;
200	<	uz = uz / u;
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	<	dAtom->setSUx(ux);
257	<	dAtom->setSUy(uy);
204	<	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
#	Line 222 \| Line 276 \| void SimSetup::makeMolecules(void){
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
228	–	// exclude_I must always be the smaller of the pair
229	–	if (exI > exJ){
230	–	tempEx = exI;
231	–	exI = exJ;
232	–	exJ = tempEx;
233	–	}
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);
241	<	#else // isn't MPI
242	<
243	<	info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
244	<	#endif //is_mpi
290	>	info[k].excludes->addPair(exI, exJ);
291		}
246	–	excludeOffset += molInfo.nBonds;
292
293		//make the bends
294		for (j = 0; j < molInfo.nBends; j++){
#	Line 293 \| Line 338 \| void SimSetup::makeMolecules(void){
338		}
339		}
340
341	<	if (!theBends[j].isGhost){
342	<	exI = theBends[j].a;
343	<	exJ = theBends[j].c;
344	<	}
345	<	else{
301	<	exI = theBends[j].a;
302	<	exJ = theBends[j].b;
303	<	}
304	<
305	<	// exclude_I must always be the smaller of the pair
306	<	if (exI > exJ){
307	<	tempEx = exI;
308	<	exI = exJ;
309	<	exJ = tempEx;
310	<	}
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		}
322	–	excludeOffset += molInfo.nBends;
376
377		for (j = 0; j < molInfo.nTorsions; j++){
378		currentTorsion = comp_stamps[stampID]->getTorsion(j);
#	Line 328 \| Line 381 \| void SimSetup::makeMolecules(void){
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
334	–	// exclude_I must always be the smaller of the pair
335	–	if (exI > exJ){
336	–	tempEx = exI;
337	–	exI = exJ;
338	–	exJ = tempEx;
339	–	}
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		}
351	–	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
362	–
478		info[k].molecules[i].initialize(molInfo);
479
480
#	Line 367 \| Line 482 \| void SimSetup::makeMolecules(void){
482		delete[] theBonds;
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
#	Line 377 \| Line 513 \| void SimSetup::makeMolecules(void){
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
#	Line 551 \| Line 695 \| void SimSetup::gatherInfo(void){
695
696
697		void SimSetup::gatherInfo(void){
698	<	int i, j, k;
698	>	int i;
699
700		ensembleCase = -1;
701		ffCase = -1;
#	Line 579 \| Line 723 \| void SimSetup::gatherInfo(void){
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);
#	Line 602 \| Line 749 \| void SimSetup::gatherInfo(void){
749		else if (!strcasecmp(ensemble, "NPTf")){
750		ensembleCase = NPTf_ENS;
751		}
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();
#	Line 637 \| Line 787 \| void SimSetup::gatherInfo(void){
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		}
#	Line 657 \| Line 807 \| void SimSetup::gatherInfo(void){
807		painCave.isFatal = 1;
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
#	Line 687 \| Line 878 \| void SimSetup::gatherInfo(void){
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
885	>	// check for the extended State init
886
887	<	double boxVector[3];
888	<	if (globals->haveBox()){
889	<	boxVector[0] = globals->getBox();
699	<	boxVector[1] = globals->getBox();
700	<	boxVector[2] = globals->getBox();
701	<
702	<	info[i].setBox(boxVector);
703	<	}
704	<	else if (globals->haveDensity()){
705	<	double vol;
706	<	vol = (double) tot_nmol / globals->getDensity();
707	<	boxVector[0] = pow(vol, (1.0 / 3.0));
708	<	boxVector[1] = boxVector[0];
709	<	boxVector[2] = boxVector[0];
710	<
711	<	info[i].setBox(boxVector);
712	<	}
713	<	else{
714	<	if (!globals->haveBoxX()){
715	<	sprintf(painCave.errMsg,
716	<	"SimSetup error, no periodic BoxX size given.\n");
717	<	painCave.isFatal = 1;
718	<	simError();
719	<	}
720	<	boxVector[0] = globals->getBoxX();
721	<
722	<	if (!globals->haveBoxY()){
723	<	sprintf(painCave.errMsg,
724	<	"SimSetup error, no periodic BoxY size given.\n");
725	<	painCave.isFatal = 1;
726	<	simError();
727	<	}
728	<	boxVector[1] = globals->getBoxY();
729	<
730	<	if (!globals->haveBoxZ()){
731	<	sprintf(painCave.errMsg,
732	<	"SimSetup error, no periodic BoxZ size given.\n");
733	<	painCave.isFatal = 1;
734	<	simError();
735	<	}
736	<	boxVector[2] = globals->getBoxZ();
737	<
738	<	info[i].setBox(boxVector);
739	<	}
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
#	Line 778 \| Line 928 \| void SimSetup::gatherInfo(void){
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		}
#	Line 813 \| Line 963 \| void SimSetup::finalInfoCheck(void){
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;
822	<	smallest = info[i].boxL[0];
823	<	if (info[i].boxL[1] <= smallest)
824	<	smallest = info[i].boxL[1];
825	<	if (info[i].boxL[2] <= smallest)
826	<	smallest = info[i].boxL[2];
827	<	theEcr = 0.5 * smallest;
971	>	theEcr = 15.0;
972		}
973		else{
974		theEcr = globals->getECR();
#	Line 832 \| Line 976 \| void SimSetup::finalInfoCheck(void){
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;
#	Line 842 \| Line 988 \| void SimSetup::finalInfoCheck(void){
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		}
#	Line 857 \| Line 1004 \| void SimSetup::finalInfoCheck(void){
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;
866	<	smallest = info[i].boxL[0];
867	<	if (info[i].boxL[1] <= smallest)
868	<	smallest = info[i].boxL[1];
869	<	if (info[i].boxL[2] <= smallest)
870	<	smallest = info[i].boxL[2];
871	<	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;
#	Line 886 \| Line 1028 \| void SimSetup::finalInfoCheck(void){
1028		else{
1029		theEst = globals->getEST();
1030		}
1031	<
1032	<	info[i].setEcr(theEcr, theEst);
1031	>
1032	>	info[i].setDefaultEcr(theEcr, theEst);
1033		}
1034		}
1035		}
894	–
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
#	Line 914 \| Line 1055 \| void SimSetup::initSystemCoords(void){
1055		if (worldRank == 0){
1056		#endif //is_mpi
1057		inName = globals->getInitialConfig();
917	–	double* tempDouble = new double[1000000];
1058		fileInit = new InitializeFromFile(inName);
1059		#ifdef IS_MPI
1060		}
#	Line 926 \| Line 1066 \| void SimSetup::initSystemCoords(void){
1066		delete fileInit;
1067		}
1068		else{
1069	<	#ifdef IS_MPI
930	<
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	<
938	<	#else
939	<
940	<	initFromBass();
941	<
942	<
943	<	#endif
1076	>
1077		}
1078
1079		#ifdef IS_MPI
#	Line 1094 \| Line 1227 \| void SimSetup::createFF(void){
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");
#	Line 1158 \| Line 1295 \| void SimSetup::calcSysValues(void){
1295		}
1296
1297		void SimSetup::calcSysValues(void){
1298	<	int i, j, k;
1298	>	int i;
1299
1300		int* molMembershipArray;
1301
#	Line 1166 \| Line 1303 \| void SimSetup::calcSysValues(void){
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
#	Line 1194 \| Line 1333 \| void SimSetup::mpiMolDivide(void){
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
#	Line 1210 \| Line 1350 \| void SimSetup::mpiMolDivide(void){
1350		local_bonds = 0;
1351		local_bends = 0;
1352		local_torsions = 0;
1353	+	local_rigid = 0;
1354		globalAtomIndex = 0;
1355
1215	–
1356		for (i = 0; i < n_components; i++){
1357		for (j = 0; j < components_nmol[i]; j++){
1358		if (mol2proc[allMol] == worldRank){
#	Line 1220 \| Line 1360 \| void SimSetup::mpiMolDivide(void){
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++){
#	Line 1233 \| Line 1374 \| void SimSetup::mpiMolDivide(void){
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();
#	Line 1257 \| Line 1399 \| void SimSetup::makeSysArrays(void){
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;
1264	–	Exclude** the_excludes;
1410
1266	–
1411		for (l = 0; l < nInfo; l++){
1412		// create the atom and short range interaction arrays
1413
#	Line 1306 \| Line 1450 \| void SimSetup::makeSysArrays(void){
1450
1451		#endif // is_mpi
1452
1453	<
1454	<	if (info[l].n_SRI){
1455	<	Exclude::createArray(info[l].n_SRI);
1312	<	the_excludes = new Exclude * [info[l].n_SRI];
1313	<	for (int ex = 0; ex < info[l].n_SRI; ex++){
1314	<	the_excludes[ex] = new Exclude(ex);
1315	<	}
1316	<	info[l].globalExcludes = new int;
1317	<	info[l].n_exclude = info[l].n_SRI;
1318	<	}
1319	<	else{
1320	<	Exclude::createArray(1);
1321	<	the_excludes = new Exclude * ;
1322	<	the_excludes[0] = new Exclude(0);
1323	<	the_excludes[0]->setPair(0, 0);
1324	<	info[l].globalExcludes = new int;
1325	<	info[l].globalExcludes[0] = 0;
1326	<	info[l].n_exclude = 0;
1327	<	}
1328	<
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;
1334	–	info[l].excludes = the_excludes;
1461
1462		the_ff->setSimInfo(info);
1463		}
#	Line 1340 \| Line 1466 \| void SimSetup::makeIntegrator(void){
1466		void SimSetup::makeIntegrator(void){
1467		int k;
1468
1469	+	NVE<RealIntegrator>* myNVE = NULL;
1470		NVT<RealIntegrator>* myNVT = 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:
#	Line 1370 \| Line 1501 \| void SimSetup::makeIntegrator(void){
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:
#	Line 1391 \| Line 1524 \| void SimSetup::makeIntegrator(void){
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		}
#	Line 1401 \| Line 1534 \| void SimSetup::makeIntegrator(void){
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		}
#	Line 1411 \| Line 1544 \| void SimSetup::makeIntegrator(void){
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:
#	Line 1432 \| Line 1567 \| void SimSetup::makeIntegrator(void){
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());
1588	+
1589		else{
1590		sprintf(painCave.errMsg,
1591		"SimSetup error: If you use an NPT\n"
1592	<	" ensemble, you must set tauBarostat.\n");
1592	>	"\tensemble, you must set tauBarostat.\n");
1593		painCave.isFatal = 1;
1594		simError();
1595		}
1596	+
1597	+	info->the_integrator = myNPTf;
1598		break;
1599
1600	+	case NPTxyz_ENS:
1601	+	if (globals->haveZconstraints()){
1602	+	setupZConstraint(info[k]);
1603	+	myNPTxyz = new ZConstraint<NPTxyz<NPT <RealIntegrator> > >(&(info[k]), the_ff);
1604	+	}
1605	+	else
1606	+	myNPTxyz = new NPTxyz<NPT <RealIntegrator> >(&(info[k]), the_ff);
1607	+
1608	+	myNPTxyz->setTargetTemp(globals->getTargetTemp());
1609	+
1610	+	if (globals->haveTargetPressure())
1611	+	myNPTxyz->setTargetPressure(globals->getTargetPressure());
1612	+	else{
1613	+	sprintf(painCave.errMsg,
1614	+	"SimSetup error: If you use a constant pressure\n"
1615	+	"\tensemble, you must set targetPressure in the BASS file.\n");
1616	+	painCave.isFatal = 1;
1617	+	simError();
1618	+	}
1619	+
1620	+	if (globals->haveTauThermostat())
1621	+	myNPTxyz->setTauThermostat(globals->getTauThermostat());
1622	+	else{
1623	+	sprintf(painCave.errMsg,
1624	+	"SimSetup error: If you use an NPT\n"
1625	+	"\tensemble, you must set tauThermostat.\n");
1626	+	painCave.isFatal = 1;
1627	+	simError();
1628	+	}
1629	+
1630	+	if (globals->haveTauBarostat())
1631	+	myNPTxyz->setTauBarostat(globals->getTauBarostat());
1632	+	else{
1633	+	sprintf(painCave.errMsg,
1634	+	"SimSetup error: If you use an NPT\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:
1644		sprintf(painCave.errMsg,
1645		"SimSetup Error. Unrecognized ensemble in case statement.\n");
#	Line 1503 \| Line 1685 \| void SimSetup::setupZConstraint(SimInfo& theInfo){
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		}
#	Line 1519 \| Line 1701 \| void SimSetup::setupZConstraint(SimInfo& theInfo){
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();
#	Line 1538 \| Line 1721 \| void SimSetup::setupZConstraint(SimInfo& theInfo){
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");
1730		}
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
#	Line 1576 \| Line 1787 \| void SimSetup::setupZConstraint(SimInfo& theInfo){
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		}
#	Line 1594 \| Line 1806 \| void SimSetup::setupZConstraint(SimInfo& theInfo){
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	+	}

Diff Legend

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

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents): Revision 780 by mmeineke, Mon Sep 22 21:23:25 2003 UTC vs. Revision 1097 by gezelter, Mon Apr 12 20:32:20 2004 UTC

Diff Legend

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 780 by mmeineke, Mon Sep 22 21:23:25 2003 UTC vs.
Revision 1097 by gezelter, Mon Apr 12 20:32:20 2004 UTC