ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/brains/ForceField.cpp
Revision: 821
Committed: Wed Dec 28 21:42:43 2005 UTC (19 years, 4 months ago) by chuckv
Original Path: trunk/src/UseTheForce/ForceField.cpp
File size: 8343 byte(s)
Log Message:
Force options passed to fortran.

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 206 #include "UseTheForce/ForceField.hpp"
51 gezelter 246 #include "utils/simError.h"
52 tim 475 #include "UseTheForce/DarkSide/atype_interface.h"
53 chuckv 821 #include "UseTheForce/DarkSide/fForceOptions_interface.h"
54 gezelter 246 namespace oopse {
55 gezelter 206
56 gezelter 507 ForceField::ForceField() {
57 gezelter 246 char* tempPath;
58     tempPath = getenv("FORCE_PARAM_PATH");
59 gezelter 206
60 gezelter 246 if (tempPath == NULL) {
61 gezelter 507 //convert a macro from compiler to a string in c++
62     STR_DEFINE(ffPath_, FRC_PATH );
63 gezelter 246 } else {
64 gezelter 507 ffPath_ = tempPath;
65 gezelter 246 }
66 gezelter 507 }
67 gezelter 206
68 tim 475
69 gezelter 507 ForceField::~ForceField() {
70 tim 475 deleteAtypes();
71 gezelter 507 }
72 tim 475
73 gezelter 507 AtomType* ForceField::getAtomType(const std::string &at) {
74 gezelter 246 std::vector<std::string> keys;
75     keys.push_back(at);
76     return atomTypeCont_.find(keys);
77 gezelter 507 }
78 gezelter 206
79 gezelter 507 BondType* ForceField::getBondType(const std::string &at1, const std::string &at2) {
80 gezelter 246 std::vector<std::string> keys;
81     keys.push_back(at1);
82     keys.push_back(at2);
83 gezelter 206
84 gezelter 246 //try exact match first
85     BondType* bondType = bondTypeCont_.find(keys);
86     if (bondType) {
87 gezelter 507 return bondType;
88 gezelter 246 } else {
89 gezelter 507 //if no exact match found, try wild card match
90     return bondTypeCont_.find(keys, wildCardAtomTypeName_);
91 gezelter 246 }
92 gezelter 206
93 gezelter 507 }
94 gezelter 206
95 gezelter 507 BendType* ForceField::getBendType(const std::string &at1, const std::string &at2,
96     const std::string &at3) {
97 gezelter 246 std::vector<std::string> keys;
98     keys.push_back(at1);
99     keys.push_back(at2);
100     keys.push_back(at3);
101 gezelter 206
102 gezelter 246 //try exact match first
103     BendType* bendType = bendTypeCont_.find(keys);
104     if (bendType) {
105 gezelter 507 return bendType;
106 gezelter 246 } else {
107 gezelter 507 //if no exact match found, try wild card match
108     return bendTypeCont_.find(keys, wildCardAtomTypeName_);
109 gezelter 246 }
110 gezelter 507 }
111 gezelter 206
112 gezelter 507 TorsionType* ForceField::getTorsionType(const std::string &at1, const std::string &at2,
113     const std::string &at3, const std::string &at4) {
114 gezelter 246 std::vector<std::string> keys;
115     keys.push_back(at1);
116     keys.push_back(at2);
117     keys.push_back(at3);
118     keys.push_back(at4);
119 gezelter 206
120 gezelter 246 TorsionType* torsionType = torsionTypeCont_.find(keys);
121     if (torsionType) {
122 gezelter 507 return torsionType;
123 gezelter 246 } else {
124 gezelter 507 //if no exact match found, try wild card match
125     return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
126 gezelter 246 }
127 gezelter 206
128 gezelter 246 return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
129 gezelter 206
130 gezelter 507 }
131 gezelter 206
132 gezelter 507 BondType* ForceField::getExactBondType(const std::string &at1, const std::string &at2){
133 gezelter 246 std::vector<std::string> keys;
134     keys.push_back(at1);
135     keys.push_back(at2);
136     return bondTypeCont_.find(keys);
137 gezelter 507 }
138 gezelter 206
139 gezelter 507 BendType* ForceField::getExactBendType(const std::string &at1, const std::string &at2,
140     const std::string &at3){
141 gezelter 246 std::vector<std::string> keys;
142     keys.push_back(at1);
143     keys.push_back(at2);
144     keys.push_back(at3);
145     return bendTypeCont_.find(keys);
146 gezelter 507 }
147 gezelter 206
148 gezelter 507 TorsionType* ForceField::getExactTorsionType(const std::string &at1, const std::string &at2,
149     const std::string &at3, const std::string &at4){
150 gezelter 246 std::vector<std::string> keys;
151     keys.push_back(at1);
152     keys.push_back(at2);
153     keys.push_back(at3);
154     keys.push_back(at4);
155     return torsionTypeCont_.find(keys);
156 gezelter 507 }
157     bool ForceField::addAtomType(const std::string &at, AtomType* atomType) {
158 gezelter 246 std::vector<std::string> keys;
159     keys.push_back(at);
160     return atomTypeCont_.add(keys, atomType);
161 gezelter 507 }
162 gezelter 206
163 gezelter 507 bool ForceField::addBondType(const std::string &at1, const std::string &at2, BondType* bondType) {
164 gezelter 246 std::vector<std::string> keys;
165     keys.push_back(at1);
166     keys.push_back(at2);
167     return bondTypeCont_.add(keys, bondType);
168 gezelter 206
169 gezelter 507 }
170 gezelter 206
171 gezelter 507 bool ForceField::addBendType(const std::string &at1, const std::string &at2,
172     const std::string &at3, BendType* bendType) {
173 gezelter 246 std::vector<std::string> keys;
174     keys.push_back(at1);
175     keys.push_back(at2);
176     keys.push_back(at3);
177     return bendTypeCont_.add(keys, bendType);
178 gezelter 507 }
179 gezelter 206
180 gezelter 507 bool ForceField::addTorsionType(const std::string &at1, const std::string &at2,
181     const std::string &at3, const std::string &at4, TorsionType* torsionType) {
182 gezelter 246 std::vector<std::string> keys;
183     keys.push_back(at1);
184     keys.push_back(at2);
185     keys.push_back(at3);
186     keys.push_back(at4);
187     return torsionTypeCont_.add(keys, torsionType);
188 gezelter 507 }
189 gezelter 206
190 gezelter 507 double ForceField::getRcutFromAtomType(AtomType* at) {
191 gezelter 246 /**@todo */
192     GenericData* data;
193     double rcut = 0.0;
194 gezelter 206
195 gezelter 246 if (at->isLennardJones()) {
196 gezelter 507 data = at->getPropertyByName("LennardJones");
197     if (data != NULL) {
198     LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
199 gezelter 206
200 gezelter 507 if (ljData != NULL) {
201     LJParam ljParam = ljData->getData();
202 gezelter 206
203 gezelter 507 //by default use 2.5*sigma as cutoff radius
204     rcut = 2.5 * ljParam.sigma;
205 gezelter 246
206 gezelter 507 } else {
207     sprintf( painCave.errMsg,
208     "Can not cast GenericData to LJParam\n");
209     painCave.severity = OOPSE_ERROR;
210     painCave.isFatal = 1;
211     simError();
212     }
213     } else {
214     sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n");
215     painCave.severity = OOPSE_ERROR;
216     painCave.isFatal = 1;
217     simError();
218     }
219 gezelter 246 }
220 gezelter 206
221 gezelter 246 return rcut;
222 gezelter 507 }
223 gezelter 206
224 gezelter 246
225 gezelter 507 ifstrstream* ForceField::openForceFieldFile(const std::string& filename) {
226 gezelter 246 std::string forceFieldFilename(filename);
227     ifstrstream* ffStream = new ifstrstream();
228    
229     //try to open the force filed file in current directory first
230     ffStream->open(forceFieldFilename.c_str());
231     if(!ffStream->is_open()){
232    
233 gezelter 507 forceFieldFilename = ffPath_ + "/" + forceFieldFilename;
234     ffStream->open( forceFieldFilename.c_str() );
235 gezelter 246
236 gezelter 507 //if current directory does not contain the force field file,
237     //try to open it in the path
238     if(!ffStream->is_open()){
239 gezelter 246
240 gezelter 507 sprintf( painCave.errMsg,
241     "Error opening the force field parameter file:\n"
242     "\t%s\n"
243     "\tHave you tried setting the FORCE_PARAM_PATH environment "
244     "variable?\n",
245     forceFieldFilename.c_str() );
246     painCave.severity = OOPSE_ERROR;
247     painCave.isFatal = 1;
248     simError();
249     }
250 gezelter 246 }
251    
252     return ffStream;
253    
254 gezelter 507 }
255 gezelter 246
256 chuckv 821 void ForceField::setFortranForceOptions(){
257     ForceOptions theseFortranOptions;
258     forceFieldOptions_.makeFortranOptions(theseFortranOptions);
259     setfForceOptions(&theseFortranOptions);
260     }
261 gezelter 246 } //end namespace oopse