ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-4/src/UseTheForce/ForceField.cpp
(Generate patch)

Comparing branches/new_design/OOPSE-4/src/UseTheForce/ForceField.cpp (file contents):
Revision 1716, Fri Nov 5 21:04:13 2004 UTC vs.
Revision 1847 by tim, Sat Dec 4 05:24:07 2004 UTC

# Line 1 | Line 1
1 + /*
2 + * Copyright (C) 2000-2004  Object Oriented Parallel Simulation Engine (OOPSE) project
3 + *
4 + * Contact: oopse@oopse.org
5 + *
6 + * This program is free software; you can redistribute it and/or
7 + * modify it under the terms of the GNU Lesser General Public License
8 + * as published by the Free Software Foundation; either version 2.1
9 + * of the License, or (at your option) any later version.
10 + * All we ask is that proper credit is given for our work, which includes
11 + * - but is not limited to - adding the above copyright notice to the beginning
12 + * of your source code files, and to any copyright notice that you may distribute
13 + * with programs based on this work.
14 + *
15 + * This program is distributed in the hope that it will be useful,
16 + * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18 + * GNU Lesser General Public License for more details.
19 + *
20 + * You should have received a copy of the GNU Lesser General Public License
21 + * along with this program; if not, write to the Free Software
22 + * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
23 + *
24 + */
25 +
26 + /**
27 +  * @file ForceField.cpp
28 +  * @author tlin
29 +  * @date 11/04/2004
30 +  * @time 22:51am
31 +  * @version 1.0
32 +  */
33 +  
34   #include "UseTheForce/ForceField.hpp"
35 + #include "utils/simError.h"
36 + namespace oopse {
37  
38 < AtomType* ForceField::getMatchingAtomType(const string &at) {
38 > ForceField::ForceField() {
39 >    char* tempPath;
40 >    tempPath = getenv("FORCE_PARAM_PATH");
41  
42 <  map<string, AtomType*>::iterator iter;
43 <  
44 <  iter = atomTypeMap.find(at);
45 <  if (iter != atomTypeMap.end()) {
46 <    return iter->second;
47 <  } else {
11 <    return NULL;
12 <  }
42 >    if (tempPath == NULL) {
43 >        //convert a macro from compiler to a string in c++
44 >        STR_DEFINE(ffPath_, FRC_PATH );
45 >    } else {
46 >        ffPath_ = tempPath;
47 >    }
48   }
49  
50 < BondType* ForceField::getMatchingBondType(const string &at1,
51 <                                          const string &at2) {
50 > AtomType* ForceField::getAtomType(const std::string &at) {
51 >    std::vector<std::string> keys;
52 >    keys.push_back(at);
53 >    return atomTypeCont_.find(keys);
54 > }
55  
56 <  map<pair<string,string>, BondType*>::iterator iter;
57 <  vector<BondType*> foundTypes;
56 > BondType* ForceField::getBondType(const std::string &at1, const std::string &at2) {
57 >    std::vector<std::string> keys;
58 >    keys.push_back(at1);
59 >    keys.push_back(at2);    
60  
61 <  iter = bondTypeMap.find(pair<at1, at2>);
62 <  if (iter != bondTypeMap.end()) {
63 <    // exact match, so just return it
64 <    return iter->second;
65 <  }
61 >    //try exact match first
62 >    BondType* bondType = bondTypeCont_.find(keys);
63 >    if (bondType) {
64 >        return bondType;
65 >    } else {
66 >        //if no exact match found, try wild card match
67 >        return bondTypeCont_.find(keys, wildCardAtomTypeName_);
68 >    }
69  
70 <  iter = bondTypeMap.find(pair<at2, at1>);
28 <  if (iter != bondTypeMap.end()) {
29 <    // exact match in reverse order, so just return it
30 <    return iter->second;
31 <  }
70 > }
71  
72 <  iter = bondTypeMap.find(pair<at1, wildCardAtomTypeName>);
73 <  if (iter != bondTypeMap.end()) {
74 <    foundTypes.push_back(iter->second);
75 <  }
72 > BendType* ForceField::getBendType(const std::string &at1, const std::string &at2,
73 >                        const std::string &at3) {
74 >    std::vector<std::string> keys;
75 >    keys.push_back(at1);
76 >    keys.push_back(at2);    
77 >    keys.push_back(at3);    
78  
79 <  iter = bondTypeMap.find(pair<at2, wildCardAtomTypeName>);
80 <  if (iter != bondTypeMap.end()) {
81 <    foundTypes.push_back(iter->second);
82 <  }
79 >    //try exact match first
80 >    BendType* bendType = bendTypeCont_.find(keys);
81 >    if (bendType) {
82 >        return bendType;
83 >    } else {
84 >        //if no exact match found, try wild card match
85 >        return bendTypeCont_.find(keys, wildCardAtomTypeName_);
86 >    }
87 > }
88  
89 <  iter = bondTypeMap.find(pair<wildCardAtomTypeName, at1>);
90 <  if (iter != bondTypeMap.end()) {
91 <    foundTypes.push_back(iter->second);
92 <  }
89 > TorsionType* ForceField::getTorsionType(const std::string &at1, const std::string &at2,
90 >                              const std::string &at3, const std::string &at4) {
91 >    std::vector<std::string> keys;
92 >    keys.push_back(at1);
93 >    keys.push_back(at2);    
94 >    keys.push_back(at3);    
95 >    keys.push_back(at4);    
96  
97 <  iter = bondTypeMap.find(pair<wildCardAtomTypeName, at2>);
98 <  if (iter != bondTypeMap.end()) {
99 <    foundTypes.push_back(iter->second);
100 <  }
101 <  
102 <  if (foundTypes.empty()) {
103 <    return NULL;
55 <  } else {
97 >    TorsionType* torsionType = torsionTypeCont_.find(keys);
98 >    if (torsionType) {
99 >        return torsionType;
100 >    } else {
101 >        //if no exact match found, try wild card match
102 >        return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
103 >    }
104      
105 +    return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
106  
107 <
107 > }
108  
109 + BondType* ForceField::getExactBondType(const std::string &at1, const std::string &at2){
110 +    std::vector<std::string> keys;
111 +    keys.push_back(at1);
112 +    keys.push_back(at2);    
113 +    return bondTypeCont_.find(keys);
114 + }
115  
116 + BendType* ForceField::getExactBendType(const std::string &at1, const std::string &at2,
117 +                            const std::string &at3){
118 +    std::vector<std::string> keys;
119 +    keys.push_back(at1);
120 +    keys.push_back(at2);    
121 +    keys.push_back(at3);    
122 +    return bendTypeCont_.find(keys);
123 + }
124  
125 <  
125 > TorsionType* ForceField::getExactTorsionType(const std::string &at1, const std::string &at2,
126 >                                  const std::string &at3, const std::string &at4){
127 >    std::vector<std::string> keys;
128 >    keys.push_back(at1);
129 >    keys.push_back(at2);    
130 >    keys.push_back(at3);    
131 >    keys.push_back(at4);  
132 >    return torsionTypeCont_.find(keys);
133 > }
134 > bool ForceField::addAtomType(const std::string &at, AtomType* atomType) {
135 >    std::vector<std::string> keys;
136 >    keys.push_back(at);
137 >    return atomTypeCont_.add(keys, atomType);
138 > }
139  
140 + bool ForceField::addBondType(const std::string &at1, const std::string &at2, BondType* bondType) {
141 +    std::vector<std::string> keys;
142 +    keys.push_back(at1);
143 +    keys.push_back(at2);    
144 +    return bondTypeCont_.add(keys, bondType);
145  
146 < BendType* ForceField::getMatchingBendType(const string &at1, const string &at2,
66 <                                          const string &at3);
67 < TorsionType* ForceField::getMatchingTorsionType(const string &at1, const string &at2,
68 <                                                const string &at3, const string &at4);
146 > }
147  
148 < double ForceField::getRcutForAtomType(AtomType* at);
148 > bool ForceField::addBendType(const std::string &at1, const std::string &at2,
149 >                        const std::string &at3, BendType* bendType) {
150 >    std::vector<std::string> keys;
151 >    keys.push_back(at1);
152 >    keys.push_back(at2);    
153 >    keys.push_back(at3);    
154 >    return bendTypeCont_.add(keys, bendType);
155 > }
156  
157 + bool ForceField::addTorsionType(const std::string &at1, const std::string &at2,
158 +                              const std::string &at3, const std::string &at4, TorsionType* torsionType) {
159 +    std::vector<std::string> keys;
160 +    keys.push_back(at1);
161 +    keys.push_back(at2);    
162 +    keys.push_back(at3);    
163 +    keys.push_back(at4);    
164 +    return torsionTypeCont_.add(keys, torsionType);
165 + }
166  
167 < vector<vector<string> > generateWildcardSequence(const vector<string> atomTypes) {
168 <  
169 <   vector<vector<string> > results;
167 > double ForceField::getRcutFromAtomType(AtomType* at) {
168 >    /**@todo */
169 >    GenericData* data;
170 >    double rcut = 0.0;
171  
172 <  
172 >    if (at->isLennardJones()) {
173 >        data = at->getPropertyByName("LennardJones");
174 >        if (data != NULL) {
175 >            LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
176  
177 +            if (ljData != NULL) {
178 +                LJParam ljParam = ljData->getData();
179  
180 <   vector<vector< string> > getAllWildcardPermutations(const vector<string> myAts) {
181 <    
182 <     int nStrings;
183 <     vector<string> oneResult;
184 <     vector<vector<string> > allResults;
180 >                //by default use 2.5*sigma as cutoff radius
181 >                rcut = 2.5 * ljParam.sigma;
182 >                
183 >            } else {
184 >                    sprintf( painCave.errMsg,
185 >                           "Can not cast GenericData to LJParam\n");
186 >                    painCave.severity = OOPSE_ERROR;
187 >                    painCave.isFatal = 1;
188 >                    simError();          
189 >            }            
190 >        } else {
191 >            sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n");
192 >            painCave.severity = OOPSE_ERROR;
193 >            painCave.isFatal = 1;
194 >            simError();          
195 >        }
196 >    }
197  
198 <     nStrings = myAts.size();
198 >    return rcut;    
199 > }
200  
201 <     if (nStrings == 1) {
202 <       oneResult.push_back(wildcardCharacter);
203 <       allResults.push_back(oneResult);
204 <       return allResults;
205 <     } else {
206 <      
207 <       for (i=0; i < nStrings; i++) {
208 <         oneResult = myAts;
209 <         replace(oneResult.begin(), oneResult.end(),
201 >
202 > ifstrstream* ForceField::openForceFieldFile(const std::string& filename) {
203 >    std::string forceFieldFilename(filename);
204 >    ifstrstream* ffStream = new ifstrstream();
205 >    
206 >    //try to open the force filed file in current directory first    
207 >    ffStream->open(forceFieldFilename.c_str());
208 >    if(!ffStream->is_open()){
209 >
210 >        forceFieldFilename = ffPath_ + "/" + forceFieldFilename;
211 >        ffStream->open( forceFieldFilename.c_str() );
212 >
213 >        //if current directory does not contain the force field file,
214 >        //try to open it in the path        
215 >        if(!ffStream->is_open()){
216 >
217 >            sprintf( painCave.errMsg,
218 >               "Error opening the force field parameter file:\n"
219 >               "\t%s\n"
220 >               "\tHave you tried setting the FORCE_PARAM_PATH environment "
221 >               "variable?\n",
222 >               forceFieldFilename.c_str() );
223 >            painCave.severity = OOPSE_ERROR;
224 >            painCave.isFatal = 1;
225 >            simError();
226 >        }
227 >    }  
228 >
229 >    return ffStream;
230 >
231 > }
232 >
233 > } //end namespace oopse

Diff Legend

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