ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/UseTheForce/ForceField.cpp
Revision: 1665
Committed: Tue Nov 22 20:38:56 2011 UTC (13 years, 5 months ago) by gezelter
File size: 21521 byte(s)
Log Message: 
updated copyright notices

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

Properties

Name Value
svn:keywords Author Id Revision Date