ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/UseTheForce/ForceField.cpp
Revision: 1282
Committed: Wed Jul 30 18:11:19 2008 UTC (16 years, 9 months ago) by gezelter
File size: 20258 byte(s)
Log Message: 
Many fixes

File Contents

# User Rev Content
1 gezelter 507 /*
2 gezelter 246 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3     *
4     * The University of Notre Dame grants you ("Licensee") a
5     * non-exclusive, royalty free, license to use, modify and
6     * redistribute this software in source and binary code form, provided
7     * that the following conditions are met:
8     *
9     * 1. Acknowledgement of the program authors must be made in any
10     * publication of scientific results based in part on use of the
11     * program. An acceptable form of acknowledgement is citation of
12     * the article in which the program was described (Matthew
13     * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14     * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15     * Parallel Simulation Engine for Molecular Dynamics,"
16     * J. Comput. Chem. 26, pp. 252-271 (2005))
17     *
18     * 2. Redistributions of source code must retain the above copyright
19     * notice, this list of conditions and the following disclaimer.
20     *
21     * 3. Redistributions in binary form must reproduce the above copyright
22     * notice, this list of conditions and the following disclaimer in the
23     * documentation and/or other materials provided with the
24     * distribution.
25     *
26     * This software is provided "AS IS," without a warranty of any
27     * kind. All express or implied conditions, representations and
28     * warranties, including any implied warranty of merchantability,
29     * fitness for a particular purpose or non-infringement, are hereby
30     * excluded. The University of Notre Dame and its licensors shall not
31     * be liable for any damages suffered by licensee as a result of
32     * using, modifying or distributing the software or its
33     * derivatives. In no event will the University of Notre Dame or its
34     * licensors be liable for any lost revenue, profit or data, or for
35     * direct, indirect, special, consequential, incidental or punitive
36     * damages, however caused and regardless of the theory of liability,
37     * arising out of the use of or inability to use software, even if the
38     * University of Notre Dame has been advised of the possibility of
39     * such damages.
40     */
41    
42 gezelter 507 /**
43     * @file ForceField.cpp
44     * @author tlin
45     * @date 11/04/2004
46     * @time 22:51am
47     * @version 1.0
48     */
49 gezelter 246
50 gezelter 1269 #include <algorithm>
51 gezelter 206 #include "UseTheForce/ForceField.hpp"
52 gezelter 246 #include "utils/simError.h"
53 gezelter 1269 #include "utils/Tuple.hpp"
54 tim 475 #include "UseTheForce/DarkSide/atype_interface.h"
55 chuckv 821 #include "UseTheForce/DarkSide/fForceOptions_interface.h"
56 gezelter 939 #include "UseTheForce/DarkSide/switcheroo_interface.h"
57 gezelter 246 namespace oopse {
58 gezelter 206
59 gezelter 507 ForceField::ForceField() {
60 gezelter 246 char* tempPath;
61     tempPath = getenv("FORCE_PARAM_PATH");
62 gezelter 206
63 gezelter 246 if (tempPath == NULL) {
64 gezelter 507 //convert a macro from compiler to a string in c++
65     STR_DEFINE(ffPath_, FRC_PATH );
66 gezelter 246 } else {
67 gezelter 507 ffPath_ = tempPath;
68 gezelter 246 }
69 gezelter 507 }
70 gezelter 206
71 tim 475
72 gezelter 507 ForceField::~ForceField() {
73 tim 475 deleteAtypes();
74 gezelter 939 deleteSwitch();
75 gezelter 507 }
76 tim 475
77 gezelter 507 AtomType* ForceField::getAtomType(const std::string &at) {
78 gezelter 246 std::vector<std::string> keys;
79     keys.push_back(at);
80     return atomTypeCont_.find(keys);
81 gezelter 507 }
82 gezelter 206
83 cpuglis 1195 BondType* ForceField::getBondType(const std::string &at1,
84     const std::string &at2) {
85 gezelter 246 std::vector<std::string> keys;
86     keys.push_back(at1);
87     keys.push_back(at2);
88 gezelter 206
89 gezelter 246 //try exact match first
90     BondType* bondType = bondTypeCont_.find(keys);
91     if (bondType) {
92 gezelter 507 return bondType;
93 gezelter 246 } else {
94 gezelter 1269 AtomType* atype1;
95     AtomType* atype2;
96     std::vector<std::string> at1key;
97     at1key.push_back(at1);
98     atype1 = atomTypeCont_.find(at1key);
99    
100     std::vector<std::string> at2key;
101     at2key.push_back(at2);
102     atype2 = atomTypeCont_.find(at2key);
103    
104     // query atom types for their chains of responsibility
105     std::vector<AtomType*> at1Chain = atype1->allYourBase();
106     std::vector<AtomType*> at2Chain = atype2->allYourBase();
107    
108     std::vector<AtomType*>::iterator i;
109     std::vector<AtomType*>::iterator j;
110    
111     int ii = 0;
112     int jj = 0;
113     int bondTypeScore;
114    
115     std::vector<std::pair<int, std::vector<std::string> > > foundBonds;
116    
117     for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
118     jj = 0;
119     for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
120    
121     bondTypeScore = ii + jj;
122    
123     std::vector<std::string> myKeys;
124     myKeys.push_back((*i)->getName());
125     myKeys.push_back((*j)->getName());
126    
127 gezelter 1277 std::cerr << "looking for " << myKeys[0] << " " << myKeys[1] << "\n";
128 gezelter 1269 BondType* bondType = bondTypeCont_.find(myKeys);
129     if (bondType) {
130     foundBonds.push_back(std::make_pair(bondTypeScore, myKeys));
131     }
132     jj++;
133     }
134     ii++;
135     }
136    
137    
138 gezelter 1277 if (foundBonds.size() > 0) {
139     // sort the foundBonds by the score:
140     std::sort(foundBonds.begin(), foundBonds.end());
141    
142     int bestScore = foundBonds[0].first;
143     std::vector<std::string> theKeys = foundBonds[0].second;
144    
145     std::cout << "best matching bond = " << theKeys[0] << "\t" << theKeys[1] << "\t(score = "<< bestScore << ")\n";
146     BondType* bestType = bondTypeCont_.find(theKeys);
147    
148     return bestType;
149     } else {
150     //if no exact match found, try wild card match
151     return bondTypeCont_.find(keys, wildCardAtomTypeName_);
152 gezelter 1269 }
153 gezelter 246 }
154 gezelter 507 }
155 gezelter 1269
156 cpuglis 1195 BendType* ForceField::getBendType(const std::string &at1,
157     const std::string &at2,
158 gezelter 507 const std::string &at3) {
159 gezelter 246 std::vector<std::string> keys;
160     keys.push_back(at1);
161     keys.push_back(at2);
162     keys.push_back(at3);
163 gezelter 206
164 gezelter 246 //try exact match first
165     BendType* bendType = bendTypeCont_.find(keys);
166     if (bendType) {
167 gezelter 507 return bendType;
168 gezelter 246 } else {
169 gezelter 1269
170     AtomType* atype1;
171     AtomType* atype2;
172     AtomType* atype3;
173     std::vector<std::string> at1key;
174     at1key.push_back(at1);
175     atype1 = atomTypeCont_.find(at1key);
176    
177     std::vector<std::string> at2key;
178     at2key.push_back(at2);
179     atype2 = atomTypeCont_.find(at2key);
180    
181     std::vector<std::string> at3key;
182     at3key.push_back(at3);
183     atype3 = atomTypeCont_.find(at3key);
184    
185     // query atom types for their chains of responsibility
186     std::vector<AtomType*> at1Chain = atype1->allYourBase();
187     std::vector<AtomType*> at2Chain = atype2->allYourBase();
188     std::vector<AtomType*> at3Chain = atype3->allYourBase();
189    
190     std::vector<AtomType*>::iterator i;
191     std::vector<AtomType*>::iterator j;
192     std::vector<AtomType*>::iterator k;
193    
194     int ii = 0;
195     int jj = 0;
196     int kk = 0;
197     int IKscore;
198    
199     std::vector<tuple3<int, int, std::vector<std::string> > > foundBends;
200    
201     for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
202     ii = 0;
203     for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
204     kk = 0;
205     for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
206    
207     IKscore = ii + kk;
208    
209     std::vector<std::string> myKeys;
210     myKeys.push_back((*i)->getName());
211     myKeys.push_back((*j)->getName());
212     myKeys.push_back((*k)->getName());
213    
214     BendType* bendType = bendTypeCont_.find(myKeys);
215     if (bendType) {
216     foundBends.push_back( make_tuple3(jj, IKscore, myKeys) );
217     }
218     kk++;
219     }
220     ii++;
221     }
222     jj++;
223     }
224    
225 gezelter 1277 if (foundBends.size() > 0) {
226     std::sort(foundBends.begin(), foundBends.end());
227     int jscore = foundBends[0].first;
228     int ikscore = foundBends[0].second;
229     std::vector<std::string> theKeys = foundBends[0].third;
230    
231     std::cout << "best matching bend = " << theKeys[0] << "\t" <<theKeys[1] << "\t" << theKeys[2] << "\t(scores = "<< jscore << "\t" << ikscore << ")\n";
232    
233     BendType* bestType = bendTypeCont_.find(theKeys);
234     return bestType;
235     } else {
236 gezelter 1269 //if no exact match found, try wild card match
237     return bendTypeCont_.find(keys, wildCardAtomTypeName_);
238     }
239 gezelter 246 }
240 gezelter 507 }
241 gezelter 206
242 cpuglis 1195 TorsionType* ForceField::getTorsionType(const std::string &at1,
243     const std::string &at2,
244     const std::string &at3,
245     const std::string &at4) {
246 gezelter 246 std::vector<std::string> keys;
247     keys.push_back(at1);
248     keys.push_back(at2);
249     keys.push_back(at3);
250     keys.push_back(at4);
251 gezelter 206
252 gezelter 1269
253     //try exact match first
254 gezelter 246 TorsionType* torsionType = torsionTypeCont_.find(keys);
255     if (torsionType) {
256 gezelter 507 return torsionType;
257 gezelter 246 } else {
258 gezelter 1269
259     AtomType* atype1;
260     AtomType* atype2;
261     AtomType* atype3;
262     AtomType* atype4;
263     std::vector<std::string> at1key;
264     at1key.push_back(at1);
265     atype1 = atomTypeCont_.find(at1key);
266    
267     std::vector<std::string> at2key;
268     at2key.push_back(at2);
269     atype2 = atomTypeCont_.find(at2key);
270    
271     std::vector<std::string> at3key;
272     at3key.push_back(at3);
273     atype3 = atomTypeCont_.find(at3key);
274    
275     std::vector<std::string> at4key;
276     at4key.push_back(at4);
277     atype4 = atomTypeCont_.find(at4key);
278    
279     // query atom types for their chains of responsibility
280     std::vector<AtomType*> at1Chain = atype1->allYourBase();
281     std::vector<AtomType*> at2Chain = atype2->allYourBase();
282     std::vector<AtomType*> at3Chain = atype3->allYourBase();
283     std::vector<AtomType*> at4Chain = atype4->allYourBase();
284    
285     std::vector<AtomType*>::iterator i;
286     std::vector<AtomType*>::iterator j;
287     std::vector<AtomType*>::iterator k;
288     std::vector<AtomType*>::iterator l;
289    
290     int ii = 0;
291     int jj = 0;
292     int kk = 0;
293     int ll = 0;
294     int ILscore;
295     int JKscore;
296    
297     std::vector<tuple3<int, int, std::vector<std::string> > > foundTorsions;
298    
299     for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
300     kk = 0;
301     for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
302     ii = 0;
303     for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
304     ll = 0;
305     for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
306    
307     ILscore = ii + ll;
308     JKscore = jj + kk;
309    
310     std::vector<std::string> myKeys;
311     myKeys.push_back((*i)->getName());
312     myKeys.push_back((*j)->getName());
313     myKeys.push_back((*k)->getName());
314     myKeys.push_back((*l)->getName());
315    
316     TorsionType* torsionType = torsionTypeCont_.find(myKeys);
317     if (torsionType) {
318     foundTorsions.push_back( make_tuple3(JKscore, ILscore, myKeys) );
319     }
320     ll++;
321     }
322     ii++;
323     }
324     kk++;
325     }
326     jj++;
327     }
328    
329 gezelter 1277 if (foundTorsions.size() > 0) {
330     std::sort(foundTorsions.begin(), foundTorsions.end());
331     int jkscore = foundTorsions[0].first;
332     int ilscore = foundTorsions[0].second;
333     std::vector<std::string> theKeys = foundTorsions[0].third;
334    
335     std::cout << "best matching torsion = " << theKeys[0] << "\t" <<theKeys[1] << "\t" << theKeys[2] << "\t" << theKeys[3] << "\t(scores = "<< jkscore << "\t" << ilscore << ")\n";
336    
337     TorsionType* bestType = torsionTypeCont_.find(theKeys);
338     return bestType;
339 gezelter 1269 } else {
340     //if no exact match found, try wild card match
341     return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
342     }
343 gezelter 246 }
344 gezelter 507 }
345 gezelter 206
346 cli2 1275 InversionType* ForceField::getInversionType(const std::string &at1,
347     const std::string &at2,
348     const std::string &at3,
349     const std::string &at4) {
350     std::vector<std::string> keys;
351     keys.push_back(at1);
352     keys.push_back(at2);
353     keys.push_back(at3);
354     keys.push_back(at4);
355    
356     //try exact match first
357     InversionType* inversionType = inversionTypeCont_.find(keys);
358     if (inversionType) {
359     return inversionType;
360     } else {
361    
362     AtomType* atype1;
363     AtomType* atype2;
364     AtomType* atype3;
365     AtomType* atype4;
366     std::vector<std::string> at1key;
367     at1key.push_back(at1);
368     atype1 = atomTypeCont_.find(at1key);
369    
370     std::vector<std::string> at2key;
371     at2key.push_back(at2);
372     atype2 = atomTypeCont_.find(at2key);
373    
374     std::vector<std::string> at3key;
375     at3key.push_back(at3);
376     atype3 = atomTypeCont_.find(at3key);
377    
378     std::vector<std::string> at4key;
379     at4key.push_back(at4);
380     atype4 = atomTypeCont_.find(at4key);
381    
382     // query atom types for their chains of responsibility
383     std::vector<AtomType*> at1Chain = atype1->allYourBase();
384     std::vector<AtomType*> at2Chain = atype2->allYourBase();
385     std::vector<AtomType*> at3Chain = atype3->allYourBase();
386     std::vector<AtomType*> at4Chain = atype4->allYourBase();
387    
388     std::vector<AtomType*>::iterator i;
389     std::vector<AtomType*>::iterator j;
390     std::vector<AtomType*>::iterator k;
391     std::vector<AtomType*>::iterator l;
392    
393     int ii = 0;
394     int jj = 0;
395     int kk = 0;
396     int ll = 0;
397     int Iscore;
398     int JKLscore;
399    
400     std::vector<tuple3<int, int, std::vector<std::string> > > foundInversions;
401    
402     for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
403     kk = 0;
404     for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
405     ii = 0;
406     for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
407     ll = 0;
408     for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
409    
410     Iscore = ii;
411     JKLscore = jj + kk + ll;
412    
413     std::vector<std::string> myKeys;
414     myKeys.push_back((*i)->getName());
415     myKeys.push_back((*j)->getName());
416     myKeys.push_back((*k)->getName());
417     myKeys.push_back((*l)->getName());
418    
419     InversionType* inversionType = inversionTypeCont_.find(myKeys);
420     if (inversionType) {
421     foundInversions.push_back( make_tuple3(Iscore, JKLscore, myKeys) );
422     }
423     ll++;
424     }
425     ii++;
426     }
427     kk++;
428     }
429     jj++;
430     }
431 gezelter 1277
432     if (foundInversions.size() > 0) {
433     std::sort(foundInversions.begin(), foundInversions.end());
434     int iscore = foundInversions[0].first;
435     int jklscore = foundInversions[0].second;
436     std::vector<std::string> theKeys = foundInversions[0].third;
437    
438     std::cout << "best matching inversion = " << theKeys[0] << "\t" <<theKeys[1] << "\t" << theKeys[2] << "\t" << theKeys[3] << "\t(scores = "<< iscore << "\t" << jklscore << ")\n";
439    
440     InversionType* bestType = inversionTypeCont_.find(theKeys);
441     return bestType;
442 cli2 1275 } else {
443     //if no exact match found, try wild card match
444     return inversionTypeCont_.find(keys, wildCardAtomTypeName_);
445     }
446     }
447     }
448    
449 chuckv 1151 NonBondedInteractionType* ForceField::getNonBondedInteractionType(const std::string &at1, const std::string &at2) {
450     std::vector<std::string> keys;
451     keys.push_back(at1);
452     keys.push_back(at2);
453 cpuglis 1195
454 chuckv 1151 //try exact match first
455     NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(keys);
456     if (nbiType) {
457     return nbiType;
458     } else {
459     //if no exact match found, try wild card match
460     return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_);
461 cpuglis 1195 }
462 chuckv 1151 }
463 cpuglis 1195
464     BondType* ForceField::getExactBondType(const std::string &at1,
465     const std::string &at2){
466 gezelter 246 std::vector<std::string> keys;
467     keys.push_back(at1);
468     keys.push_back(at2);
469     return bondTypeCont_.find(keys);
470 gezelter 507 }
471 cpuglis 1195
472     BendType* ForceField::getExactBendType(const std::string &at1,
473     const std::string &at2,
474 gezelter 507 const std::string &at3){
475 gezelter 246 std::vector<std::string> keys;
476     keys.push_back(at1);
477     keys.push_back(at2);
478     keys.push_back(at3);
479     return bendTypeCont_.find(keys);
480 gezelter 507 }
481 cpuglis 1195
482     TorsionType* ForceField::getExactTorsionType(const std::string &at1,
483     const std::string &at2,
484     const std::string &at3,
485     const std::string &at4){
486 gezelter 246 std::vector<std::string> keys;
487     keys.push_back(at1);
488     keys.push_back(at2);
489     keys.push_back(at3);
490     keys.push_back(at4);
491     return torsionTypeCont_.find(keys);
492 gezelter 507 }
493 cli2 1275
494     InversionType* ForceField::getExactInversionType(const std::string &at1,
495     const std::string &at2,
496     const std::string &at3,
497     const std::string &at4){
498     std::vector<std::string> keys;
499     keys.push_back(at1);
500     keys.push_back(at2);
501     keys.push_back(at3);
502     keys.push_back(at4);
503     return inversionTypeCont_.find(keys);
504     }
505    
506 chuckv 1151 NonBondedInteractionType* ForceField::getExactNonBondedInteractionType(const std::string &at1, const std::string &at2){
507     std::vector<std::string> keys;
508     keys.push_back(at1);
509     keys.push_back(at2);
510     return nonBondedInteractionTypeCont_.find(keys);
511     }
512 cli2 1275
513 chuckv 1151
514 gezelter 507 bool ForceField::addAtomType(const std::string &at, AtomType* atomType) {
515 gezelter 246 std::vector<std::string> keys;
516     keys.push_back(at);
517     return atomTypeCont_.add(keys, atomType);
518 gezelter 507 }
519 gezelter 206
520 gezelter 1282 bool ForceField::replaceAtomType(const std::string &at, AtomType* atomType) {
521     std::vector<std::string> keys;
522     keys.push_back(at);
523     return atomTypeCont_.replace(keys, atomType);
524     }
525    
526 cpuglis 1195 bool ForceField::addBondType(const std::string &at1, const std::string &at2,
527     BondType* bondType) {
528 gezelter 246 std::vector<std::string> keys;
529     keys.push_back(at1);
530     keys.push_back(at2);
531 cpuglis 1195 return bondTypeCont_.add(keys, bondType);
532 gezelter 507 }
533 cpuglis 1195
534 gezelter 507 bool ForceField::addBendType(const std::string &at1, const std::string &at2,
535     const std::string &at3, BendType* bendType) {
536 gezelter 246 std::vector<std::string> keys;
537     keys.push_back(at1);
538     keys.push_back(at2);
539     keys.push_back(at3);
540     return bendTypeCont_.add(keys, bendType);
541 gezelter 507 }
542 cpuglis 1195
543     bool ForceField::addTorsionType(const std::string &at1,
544     const std::string &at2,
545     const std::string &at3,
546     const std::string &at4,
547     TorsionType* torsionType) {
548 gezelter 246 std::vector<std::string> keys;
549     keys.push_back(at1);
550     keys.push_back(at2);
551     keys.push_back(at3);
552     keys.push_back(at4);
553     return torsionTypeCont_.add(keys, torsionType);
554 gezelter 507 }
555 gezelter 206
556 cli2 1275 bool ForceField::addInversionType(const std::string &at1,
557     const std::string &at2,
558     const std::string &at3,
559     const std::string &at4,
560     InversionType* inversionType) {
561     std::vector<std::string> keys;
562     keys.push_back(at1);
563     keys.push_back(at2);
564     keys.push_back(at3);
565     keys.push_back(at4);
566     return inversionTypeCont_.add(keys, inversionType);
567     }
568    
569 cpuglis 1195 bool ForceField::addNonBondedInteractionType(const std::string &at1,
570     const std::string &at2,
571     NonBondedInteractionType* nbiType) {
572 chuckv 1151 std::vector<std::string> keys;
573     keys.push_back(at1);
574     keys.push_back(at2);
575     return nonBondedInteractionTypeCont_.add(keys, nbiType);
576     }
577 cpuglis 1195
578 tim 963 RealType ForceField::getRcutFromAtomType(AtomType* at) {
579 gezelter 246 /**@todo */
580     GenericData* data;
581 tim 963 RealType rcut = 0.0;
582 cpuglis 1195
583 gezelter 246 if (at->isLennardJones()) {
584 gezelter 507 data = at->getPropertyByName("LennardJones");
585     if (data != NULL) {
586     LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
587 cpuglis 1195
588 gezelter 507 if (ljData != NULL) {
589     LJParam ljParam = ljData->getData();
590 cpuglis 1195
591 gezelter 507 //by default use 2.5*sigma as cutoff radius
592     rcut = 2.5 * ljParam.sigma;
593 cpuglis 1195
594 gezelter 507 } else {
595     sprintf( painCave.errMsg,
596     "Can not cast GenericData to LJParam\n");
597     painCave.severity = OOPSE_ERROR;
598     painCave.isFatal = 1;
599     simError();
600     }
601     } else {
602     sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n");
603     painCave.severity = OOPSE_ERROR;
604     painCave.isFatal = 1;
605     simError();
606     }
607 gezelter 246 }
608     return rcut;
609 gezelter 507 }
610 cpuglis 1195
611 gezelter 206
612 gezelter 507 ifstrstream* ForceField::openForceFieldFile(const std::string& filename) {
613 gezelter 246 std::string forceFieldFilename(filename);
614     ifstrstream* ffStream = new ifstrstream();
615    
616     //try to open the force filed file in current directory first
617     ffStream->open(forceFieldFilename.c_str());
618     if(!ffStream->is_open()){
619    
620 gezelter 507 forceFieldFilename = ffPath_ + "/" + forceFieldFilename;
621     ffStream->open( forceFieldFilename.c_str() );
622 gezelter 246
623 gezelter 507 //if current directory does not contain the force field file,
624     //try to open it in the path
625     if(!ffStream->is_open()){
626 gezelter 246
627 gezelter 507 sprintf( painCave.errMsg,
628     "Error opening the force field parameter file:\n"
629     "\t%s\n"
630     "\tHave you tried setting the FORCE_PARAM_PATH environment "
631     "variable?\n",
632     forceFieldFilename.c_str() );
633     painCave.severity = OOPSE_ERROR;
634     painCave.isFatal = 1;
635     simError();
636     }
637 gezelter 246 }
638     return ffStream;
639 gezelter 507 }
640 gezelter 246
641 chuckv 821 void ForceField::setFortranForceOptions(){
642     ForceOptions theseFortranOptions;
643     forceFieldOptions_.makeFortranOptions(theseFortranOptions);
644     setfForceOptions(&theseFortranOptions);
645     }
646 gezelter 246 } //end namespace oopse