ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/OpenMD/trunk/src/brains/ForceField.cpp

Comparing:
trunk/src/UseTheForce/ForceField.cpp (file contents), Revision 206 by gezelter, Thu Nov 4 20:51:23 2004 UTC vs.
trunk/src/brains/ForceField.cpp (file contents), Revision 1880 by gezelter, Mon Jun 17 18:28:30 2013 UTC

#	Line 1 | Line 1
1	<	#include "UseTheForce/ForceField.hpp"
1	>	/*
2	>	* 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. Redistributions of source code must retain the above copyright
10	>	*    notice, this list of conditions and the following disclaimer.
11	>	*
12	>	* 2. Redistributions in binary form must reproduce the above copyright
13	>	*    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	>	*
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, 234107 (2008).
39	>	* [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40	>	* [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41	>	*/
42	>
43	>	/**
44	>	* @file ForceField.cpp
45	>	* @author tlin
46	>	* @date 11/04/2004
47	>	* @version 1.0
48	>	*/
49	>
50	>	#include <algorithm>
51	>	#include "brains/ForceField.hpp"
52	>	#include "utils/simError.h"
53
54	<	AtomType* ForceField::getMatchingAtomType(const string &at) {
54	>	#include "io/OptionSectionParser.hpp"
55	>	#include "io/BaseAtomTypesSectionParser.hpp"
56	>	#include "io/DirectionalAtomTypesSectionParser.hpp"
57	>	#include "io/AtomTypesSectionParser.hpp"
58	>	#include "io/BendTypesSectionParser.hpp"
59	>	#include "io/BondTypesSectionParser.hpp"
60	>	#include "io/ChargeAtomTypesSectionParser.hpp"
61	>	#include "io/EAMAtomTypesSectionParser.hpp"
62	>	#include "io/FluctuatingChargeAtomTypesSectionParser.hpp"
63	>	#include "io/GayBerneAtomTypesSectionParser.hpp"
64	>	#include "io/InversionTypesSectionParser.hpp"
65	>	#include "io/LennardJonesAtomTypesSectionParser.hpp"
66	>	#include "io/MultipoleAtomTypesSectionParser.hpp"
67	>	#include "io/NonBondedInteractionsSectionParser.hpp"
68	>	#include "io/PolarizableAtomTypesSectionParser.hpp"
69	>	#include "io/SCAtomTypesSectionParser.hpp"
70	>	#include "io/ShapeAtomTypesSectionParser.hpp"
71	>	#include "io/StickyAtomTypesSectionParser.hpp"
72	>	#include "io/StickyPowerAtomTypesSectionParser.hpp"
73	>	#include "io/TorsionTypesSectionParser.hpp"
74
75	<	map<string, AtomType*>::iterator iter;
76	<
77	<	iter = atomTypeMap.find(at);
78	<	if (iter != atomTypeMap.end()) {
79	<	return iter->second;
10	<	} else {
11	<	return NULL;
12	<	}
13	<	}
75	>	#include "types/LennardJonesAdapter.hpp"
76	>	#include "types/EAMAdapter.hpp"
77	>	#include "types/SuttonChenAdapter.hpp"
78	>	#include "types/GayBerneAdapter.hpp"
79	>	#include "types/StickyAdapter.hpp"
80
81	<	BondType* ForceField::getMatchingBondType(const string &at1,
16	<	const string &at2) {
81	>	namespace OpenMD {
82
83	<	map<pair<string,string>, BondType*>::iterator iter;
19	<	vector<BondType*> foundTypes;
83	>	ForceField::ForceField(std::string ffName) {
84
85	<	iter = bondTypeMap.find(pair<at1, at2>);
86	<	if (iter != bondTypeMap.end()) {
87	<	// exact match, so just return it
88	<	return iter->second;
89	<	}
85	>	char* tempPath;
86	>	tempPath = getenv("FORCE_PARAM_PATH");
87	>
88	>	if (tempPath == NULL) {
89	>	//convert a macro from compiler to a string in c++
90	>	STR_DEFINE(ffPath_, FRC_PATH );
91	>	} else {
92	>	ffPath_ = tempPath;
93	>	}
94
95	<	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	<	}
95	>	setForceFieldFileName(ffName + ".frc");
96
97	<	iter = bondTypeMap.find(pair<at1, wildCardAtomTypeName>);
98	<	if (iter != bondTypeMap.end()) {
99	<	foundTypes.push_back(iter->second);
97	>	/**
98	>	* The order of adding section parsers is important.
99	>	*
100	>	* OptionSectionParser must come first to set options for other
101	>	* parsers
102	>	*
103	>	* DirectionalAtomTypesSectionParser should be added before
104	>	* AtomTypesSectionParser, and these two section parsers will
105	>	* actually create "real" AtomTypes (AtomTypesSectionParser will
106	>	* create AtomType and DirectionalAtomTypesSectionParser will
107	>	* create DirectionalAtomType, which is a subclass of AtomType and
108	>	* should come first).
109	>	*
110	>	* Other AtomTypes Section Parsers will not create the "real"
111	>	* AtomType, they only add and set some attributes of the AtomType
112	>	* (via the Adapters). Thus ordering of these is not important.
113	>	* AtomTypesSectionParser should be added before other atom type
114	>	*
115	>	* The order of BondTypesSectionParser, BendTypesSectionParser and
116	>	* TorsionTypesSectionParser, etc. are not important.
117	>	*/
118	>
119	>	spMan_.push_back(new OptionSectionParser(forceFieldOptions_));
120	>	spMan_.push_back(new BaseAtomTypesSectionParser());
121	>	spMan_.push_back(new DirectionalAtomTypesSectionParser(forceFieldOptions_));
122	>	spMan_.push_back(new AtomTypesSectionParser());
123	>
124	>	spMan_.push_back(new LennardJonesAtomTypesSectionParser(forceFieldOptions_));
125	>	spMan_.push_back(new ChargeAtomTypesSectionParser(forceFieldOptions_));
126	>	spMan_.push_back(new MultipoleAtomTypesSectionParser(forceFieldOptions_));
127	>	spMan_.push_back(new FluctuatingChargeAtomTypesSectionParser(forceFieldOptions_));
128	>	spMan_.push_back(new PolarizableAtomTypesSectionParser(forceFieldOptions_));
129	>	spMan_.push_back(new GayBerneAtomTypesSectionParser(forceFieldOptions_));
130	>	spMan_.push_back(new EAMAtomTypesSectionParser(forceFieldOptions_));
131	>	spMan_.push_back(new SCAtomTypesSectionParser(forceFieldOptions_));
132	>	spMan_.push_back(new ShapeAtomTypesSectionParser(forceFieldOptions_));
133	>	spMan_.push_back(new StickyAtomTypesSectionParser(forceFieldOptions_));
134	>	spMan_.push_back(new StickyPowerAtomTypesSectionParser(forceFieldOptions_));
135	>
136	>	spMan_.push_back(new BondTypesSectionParser(forceFieldOptions_));
137	>	spMan_.push_back(new BendTypesSectionParser(forceFieldOptions_));
138	>	spMan_.push_back(new TorsionTypesSectionParser(forceFieldOptions_));
139	>	spMan_.push_back(new InversionTypesSectionParser(forceFieldOptions_));
140	>
141	>	spMan_.push_back(new NonBondedInteractionsSectionParser(forceFieldOptions_));
142		}
143
144	<	iter = bondTypeMap.find(pair<at2, wildCardAtomTypeName>);
145	<	if (iter != bondTypeMap.end()) {
146	<	foundTypes.push_back(iter->second);
144	>	void ForceField::parse(const std::string& filename) {
145	>	ifstrstream* ffStream;
146	>
147	>	ffStream = openForceFieldFile(filename);
148	>
149	>	spMan_.parse(*ffStream, *this);
150	>
151	>	ForceField::AtomTypeContainer::MapTypeIterator i;
152	>	AtomType* at;
153	>
154	>	for (at = atomTypeCont_.beginType(i); at != NULL;
155	>	at = atomTypeCont_.nextType(i)) {
156	>
157	>	// useBase sets the responsibilities, and these have to be done
158	>	// after the atomTypes and Base types have all been scanned:
159	>
160	>	std::vector<AtomType*> ayb = at->allYourBase();
161	>	if (ayb.size() > 1) {
162	>	for (int j = ayb.size()-1; j > 0; j--) {
163	>
164	>	ayb[j-1]->useBase(ayb[j]);
165	>
166	>	}
167	>	}
168	>	}
169	>
170	>	delete ffStream;
171		}
172
173	<	iter = bondTypeMap.find(pair<wildCardAtomTypeName, at1>);
174	<	if (iter != bondTypeMap.end()) {
175	<	foundTypes.push_back(iter->second);
173	>	/**
174	>	* getAtomType by string
175	>	*
176	>	* finds the requested atom type in this force field using the string
177	>	* name of the atom type.
178	>	*/
179	>	AtomType* ForceField::getAtomType(const std::string &at) {
180	>	std::vector<std::string> keys;
181	>	keys.push_back(at);
182	>	return atomTypeCont_.find(keys);
183		}
184
185	<	iter = bondTypeMap.find(pair<wildCardAtomTypeName, at2>);
186	<	if (iter != bondTypeMap.end()) {
187	<	foundTypes.push_back(iter->second);
185	>	/**
186	>	* getAtomType by ident
187	>	*
188	>	* finds the requested atom type in this force field using the
189	>	* integer ident instead of the string name of the atom type.
190	>	*/
191	>	AtomType* ForceField::getAtomType(int ident) {
192	>	std::string at = atypeIdentToName.find(ident)->second;
193	>	return getAtomType(at);
194		}
195	+
196	+	BondType* ForceField::getBondType(const std::string &at1,
197	+	const std::string &at2) {
198	+	std::vector<std::string> keys;
199	+	keys.push_back(at1);
200	+	keys.push_back(at2);
201	+
202	+	//try exact match first
203	+	BondType* bondType = bondTypeCont_.find(keys);
204	+	if (bondType) {
205	+	return bondType;
206	+	} else {
207	+	AtomType* atype1;
208	+	AtomType* atype2;
209	+	std::vector<std::string> at1key;
210	+	at1key.push_back(at1);
211	+	atype1 = atomTypeCont_.find(at1key);
212
213	<	if (foundTypes.empty()) {
214	<	return NULL;
215	<	} else {
56	<
213	>	std::vector<std::string> at2key;
214	>	at2key.push_back(at2);
215	>	atype2 = atomTypeCont_.find(at2key);
216
217	<
217	>	// query atom types for their chains of responsibility
218	>	std::vector<AtomType*> at1Chain = atype1->allYourBase();
219	>	std::vector<AtomType*> at2Chain = atype2->allYourBase();
220
221	+	std::vector<AtomType*>::iterator i;
222	+	std::vector<AtomType*>::iterator j;
223
224	+	int ii = 0;
225	+	int jj = 0;
226	+	int bondTypeScore;
227
228	+	std::vector<std::pair<int, std::vector<std::string> > > foundBonds;
229	+
230	+	for (i = at1Chain.begin(); i != at1Chain.end(); ++i) {
231	+	jj = 0;
232	+	for (j = at2Chain.begin(); j != at2Chain.end(); ++j) {
233	+
234	+	bondTypeScore = ii + jj;
235	+
236	+	std::vector<std::string> myKeys;
237	+	myKeys.push_back((*i)->getName());
238	+	myKeys.push_back((*j)->getName());
239	+
240	+	BondType* bondType = bondTypeCont_.find(myKeys);
241	+	if (bondType) {
242	+	foundBonds.push_back(std::make_pair(bondTypeScore, myKeys));
243	+	}
244	+	jj++;
245	+	}
246	+	ii++;
247	+	}
248	+
249	+
250	+	if (!foundBonds.empty()) {
251	+	// sort the foundBonds by the score:
252	+	std::sort(foundBonds.begin(), foundBonds.end());
253	+
254	+	std::vector<std::string> theKeys = foundBonds[0].second;
255	+
256	+	BondType* bestType = bondTypeCont_.find(theKeys);
257	+
258	+	return bestType;
259	+	} else {
260	+	//if no exact match found, try wild card match
261	+	return bondTypeCont_.find(keys, wildCardAtomTypeName_);
262	+	}
263	+	}
264	+	}
265
266	+	BendType* ForceField::getBendType(const std::string &at1,
267	+	const std::string &at2,
268	+	const std::string &at3) {
269	+	std::vector<std::string> keys;
270	+	keys.push_back(at1);
271	+	keys.push_back(at2);
272	+	keys.push_back(at3);
273
274	+	//try exact match first
275	+	BendType* bendType = bendTypeCont_.find(keys);
276	+	if (bendType) {
277	+	return bendType;
278	+	} else {
279
280	<	BendType* ForceField::getMatchingBendType(const string &at1, const string &at2,
281	<	const string &at3);
282	<	TorsionType* ForceField::getMatchingTorsionType(const string &at1, const string &at2,
283	<	const string &at3, const string &at4);
280	>	AtomType* atype1;
281	>	AtomType* atype2;
282	>	AtomType* atype3;
283	>	std::vector<std::string> at1key;
284	>	at1key.push_back(at1);
285	>	atype1 = atomTypeCont_.find(at1key);
286	>
287	>	std::vector<std::string> at2key;
288	>	at2key.push_back(at2);
289	>	atype2 = atomTypeCont_.find(at2key);
290
291	<	double ForceField::getRcutForAtomType(AtomType* at);
291	>	std::vector<std::string> at3key;
292	>	at3key.push_back(at3);
293	>	atype3 = atomTypeCont_.find(at3key);
294
295	+	// query atom types for their chains of responsibility
296	+	std::vector<AtomType*> at1Chain = atype1->allYourBase();
297	+	std::vector<AtomType*> at2Chain = atype2->allYourBase();
298	+	std::vector<AtomType*> at3Chain = atype3->allYourBase();
299
300	<	vector<vector<string> > generateWildcardSequence(const vector<string> atomTypes) {
301	<
302	<	vector<vector<string> > results;
300	>	std::vector<AtomType*>::iterator i;
301	>	std::vector<AtomType*>::iterator j;
302	>	std::vector<AtomType*>::iterator k;
303
304	<
304	>	int ii = 0;
305	>	int jj = 0;
306	>	int kk = 0;
307	>	int IKscore;
308
309	+	std::vector<tuple3<int, int, std::vector<std::string> > > foundBends;
310
311	<	vector<vector< string> > getAllWildcardPermutations(const vector<string> myAts) {
312	<
313	<	int nStrings;
314	<	vector<string> oneResult;
315	<	vector<vector<string> > allResults;
311	>	for (j = at2Chain.begin(); j != at2Chain.end(); ++j) {
312	>	ii = 0;
313	>	for (i = at1Chain.begin(); i != at1Chain.end(); ++i) {
314	>	kk = 0;
315	>	for (k = at3Chain.begin(); k != at3Chain.end(); ++k) {
316	>
317	>	IKscore = ii + kk;
318
319	<	nStrings = myAts.size();
319	>	std::vector<std::string> myKeys;
320	>	myKeys.push_back((*i)->getName());
321	>	myKeys.push_back((*j)->getName());
322	>	myKeys.push_back((*k)->getName());
323
324	<	if (nStrings == 1) {
325	<	oneResult.push_back(wildcardCharacter);
326	<	allResults.push_back(oneResult);
327	<	return allResults;
328	<	} else {
329	<
330	<	for (i=0; i < nStrings; i++) {
331	<	oneResult = myAts;
332	<	replace(oneResult.begin(), oneResult.end(),
324	>	BendType* bendType = bendTypeCont_.find(myKeys);
325	>	if (bendType) {
326	>	foundBends.push_back( make_tuple3(jj, IKscore, myKeys) );
327	>	}
328	>	kk++;
329	>	}
330	>	ii++;
331	>	}
332	>	jj++;
333	>	}
334	>
335	>	if (!foundBends.empty()) {
336	>	std::sort(foundBends.begin(), foundBends.end());
337	>	std::vector<std::string> theKeys = foundBends[0].third;
338	>
339	>	BendType* bestType = bendTypeCont_.find(theKeys);
340	>	return bestType;
341	>	} else {
342	>	//if no exact match found, try wild card match
343	>	return bendTypeCont_.find(keys, wildCardAtomTypeName_);
344	>	}
345	>	}
346	>	}
347	>
348	>	TorsionType* ForceField::getTorsionType(const std::string &at1,
349	>	const std::string &at2,
350	>	const std::string &at3,
351	>	const std::string &at4) {
352	>	std::vector<std::string> keys;
353	>	keys.push_back(at1);
354	>	keys.push_back(at2);
355	>	keys.push_back(at3);
356	>	keys.push_back(at4);
357	>
358	>
359	>	//try exact match first
360	>	TorsionType* torsionType = torsionTypeCont_.find(keys);
361	>	if (torsionType) {
362	>	return torsionType;
363	>	} else {
364	>
365	>	AtomType* atype1;
366	>	AtomType* atype2;
367	>	AtomType* atype3;
368	>	AtomType* atype4;
369	>	std::vector<std::string> at1key;
370	>	at1key.push_back(at1);
371	>	atype1 = atomTypeCont_.find(at1key);
372	>
373	>	std::vector<std::string> at2key;
374	>	at2key.push_back(at2);
375	>	atype2 = atomTypeCont_.find(at2key);
376	>
377	>	std::vector<std::string> at3key;
378	>	at3key.push_back(at3);
379	>	atype3 = atomTypeCont_.find(at3key);
380	>
381	>	std::vector<std::string> at4key;
382	>	at4key.push_back(at4);
383	>	atype4 = atomTypeCont_.find(at4key);
384	>
385	>	// query atom types for their chains of responsibility
386	>	std::vector<AtomType*> at1Chain = atype1->allYourBase();
387	>	std::vector<AtomType*> at2Chain = atype2->allYourBase();
388	>	std::vector<AtomType*> at3Chain = atype3->allYourBase();
389	>	std::vector<AtomType*> at4Chain = atype4->allYourBase();
390	>
391	>	std::vector<AtomType*>::iterator i;
392	>	std::vector<AtomType*>::iterator j;
393	>	std::vector<AtomType*>::iterator k;
394	>	std::vector<AtomType*>::iterator l;
395	>
396	>	int ii = 0;
397	>	int jj = 0;
398	>	int kk = 0;
399	>	int ll = 0;
400	>	int ILscore;
401	>	int JKscore;
402	>
403	>	std::vector<tuple3<int, int, std::vector<std::string> > > foundTorsions;
404	>
405	>	for (j = at2Chain.begin(); j != at2Chain.end(); ++j) {
406	>	kk = 0;
407	>	for (k = at3Chain.begin(); k != at3Chain.end(); ++k) {
408	>	ii = 0;
409	>	for (i = at1Chain.begin(); i != at1Chain.end(); ++i) {
410	>	ll = 0;
411	>	for (l = at4Chain.begin(); l != at4Chain.end(); ++l) {
412	>
413	>	ILscore = ii + ll;
414	>	JKscore = jj + kk;
415	>
416	>	std::vector<std::string> myKeys;
417	>	myKeys.push_back((*i)->getName());
418	>	myKeys.push_back((*j)->getName());
419	>	myKeys.push_back((*k)->getName());
420	>	myKeys.push_back((*l)->getName());
421	>
422	>	TorsionType* torsionType = torsionTypeCont_.find(myKeys);
423	>	if (torsionType) {
424	>	foundTorsions.push_back( make_tuple3(JKscore, ILscore, myKeys) );
425	>	}
426	>	ll++;
427	>	}
428	>	ii++;
429	>	}
430	>	kk++;
431	>	}
432	>	jj++;
433	>	}
434	>
435	>	if (!foundTorsions.empty()) {
436	>	std::sort(foundTorsions.begin(), foundTorsions.end());
437	>	std::vector<std::string> theKeys = foundTorsions[0].third;
438	>
439	>	TorsionType* bestType = torsionTypeCont_.find(theKeys);
440	>	return bestType;
441	>	} else {
442	>	//if no exact match found, try wild card match
443	>	return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
444	>	}
445	>	}
446	>	}
447	>
448	>	InversionType* ForceField::getInversionType(const std::string &at1,
449	>	const std::string &at2,
450	>	const std::string &at3,
451	>	const std::string &at4) {
452	>	std::vector<std::string> keys;
453	>	keys.push_back(at1);
454	>	keys.push_back(at2);
455	>	keys.push_back(at3);
456	>	keys.push_back(at4);
457	>
458	>	//try exact match first
459	>	InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(keys);
460	>	if (inversionType) {
461	>	return inversionType;
462	>	} else {
463	>
464	>	AtomType* atype1;
465	>	AtomType* atype2;
466	>	AtomType* atype3;
467	>	AtomType* atype4;
468	>	std::vector<std::string> at1key;
469	>	at1key.push_back(at1);
470	>	atype1 = atomTypeCont_.find(at1key);
471	>
472	>	std::vector<std::string> at2key;
473	>	at2key.push_back(at2);
474	>	atype2 = atomTypeCont_.find(at2key);
475	>
476	>	std::vector<std::string> at3key;
477	>	at3key.push_back(at3);
478	>	atype3 = atomTypeCont_.find(at3key);
479	>
480	>	std::vector<std::string> at4key;
481	>	at4key.push_back(at4);
482	>	atype4 = atomTypeCont_.find(at4key);
483	>
484	>	// query atom types for their chains of responsibility
485	>	std::vector<AtomType*> at1Chain = atype1->allYourBase();
486	>	std::vector<AtomType*> at2Chain = atype2->allYourBase();
487	>	std::vector<AtomType*> at3Chain = atype3->allYourBase();
488	>	std::vector<AtomType*> at4Chain = atype4->allYourBase();
489	>
490	>	std::vector<AtomType*>::iterator i;
491	>	std::vector<AtomType*>::iterator j;
492	>	std::vector<AtomType*>::iterator k;
493	>	std::vector<AtomType*>::iterator l;
494	>
495	>	int ii = 0;
496	>	int jj = 0;
497	>	int kk = 0;
498	>	int ll = 0;
499	>	int Iscore;
500	>	int JKLscore;
501	>
502	>	std::vector<tuple3<int, int, std::vector<std::string> > > foundInversions;
503	>
504	>	for (j = at2Chain.begin(); j != at2Chain.end(); ++j) {
505	>	kk = 0;
506	>	for (k = at3Chain.begin(); k != at3Chain.end(); ++k) {
507	>	ii = 0;
508	>	for (i = at1Chain.begin(); i != at1Chain.end(); ++i) {
509	>	ll = 0;
510	>	for (l = at4Chain.begin(); l != at4Chain.end(); ++l) {
511	>
512	>	Iscore = ii;
513	>	JKLscore = jj + kk + ll;
514	>
515	>	std::vector<std::string> myKeys;
516	>	myKeys.push_back((*i)->getName());
517	>	myKeys.push_back((*j)->getName());
518	>	myKeys.push_back((*k)->getName());
519	>	myKeys.push_back((*l)->getName());
520	>
521	>	InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(myKeys);
522	>	if (inversionType) {
523	>	foundInversions.push_back( make_tuple3(Iscore, JKLscore, myKeys) );
524	>	}
525	>	ll++;
526	>	}
527	>	ii++;
528	>	}
529	>	kk++;
530	>	}
531	>	jj++;
532	>	}
533	>
534	>	if (!foundInversions.empty()) {
535	>	std::sort(foundInversions.begin(), foundInversions.end());
536	>	std::vector<std::string> theKeys = foundInversions[0].third;
537	>
538	>	InversionType* bestType = inversionTypeCont_.permutedFindSkippingFirstElement(theKeys);
539	>	return bestType;
540	>	} else {
541	>	//if no exact match found, try wild card match
542	>	return inversionTypeCont_.find(keys, wildCardAtomTypeName_);
543	>	}
544	>	}
545	>	}
546	>
547	>	NonBondedInteractionType* ForceField::getNonBondedInteractionType(const std::string &at1, const std::string &at2) {
548	>
549	>	std::vector<std::string> keys;
550	>	keys.push_back(at1);
551	>	keys.push_back(at2);
552	>
553	>	//try exact match first
554	>	NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(keys);
555	>	if (nbiType) {
556	>	return nbiType;
557	>	} else {
558	>	AtomType* atype1;
559	>	AtomType* atype2;
560	>	std::vector<std::string> at1key;
561	>	at1key.push_back(at1);
562	>	atype1 = atomTypeCont_.find(at1key);
563	>
564	>	std::vector<std::string> at2key;
565	>	at2key.push_back(at2);
566	>	atype2 = atomTypeCont_.find(at2key);
567	>
568	>	// query atom types for their chains of responsibility
569	>	std::vector<AtomType*> at1Chain = atype1->allYourBase();
570	>	std::vector<AtomType*> at2Chain = atype2->allYourBase();
571	>
572	>	std::vector<AtomType*>::iterator i;
573	>	std::vector<AtomType*>::iterator j;
574	>
575	>	int ii = 0;
576	>	int jj = 0;
577	>	int nbiTypeScore;
578	>
579	>	std::vector<std::pair<int, std::vector<std::string> > > foundNBI;
580	>
581	>	for (i = at1Chain.begin(); i != at1Chain.end(); ++i) {
582	>	jj = 0;
583	>	for (j = at2Chain.begin(); j != at2Chain.end(); ++j) {
584	>
585	>	nbiTypeScore = ii + jj;
586	>
587	>	std::vector<std::string> myKeys;
588	>	myKeys.push_back((*i)->getName());
589	>	myKeys.push_back((*j)->getName());
590	>
591	>	NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(myKeys);
592	>	if (nbiType) {
593	>	foundNBI.push_back(std::make_pair(nbiTypeScore, myKeys));
594	>	}
595	>	jj++;
596	>	}
597	>	ii++;
598	>	}
599	>
600	>
601	>	if (!foundNBI.empty()) {
602	>	// sort the foundNBI by the score:
603	>	std::sort(foundNBI.begin(), foundNBI.end());
604	>	std::vector<std::string> theKeys = foundNBI[0].second;
605	>
606	>	NonBondedInteractionType* bestType = nonBondedInteractionTypeCont_.find(theKeys);
607	>	return bestType;
608	>	} else {
609	>	//if no exact match found, try wild card match
610	>	return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_);
611	>	}
612	>	}
613	>	}
614	>
615	>	BondType* ForceField::getExactBondType(const std::string &at1,
616	>	const std::string &at2){
617	>	std::vector<std::string> keys;
618	>	keys.push_back(at1);
619	>	keys.push_back(at2);
620	>	return bondTypeCont_.find(keys);
621	>	}
622	>
623	>	BendType* ForceField::getExactBendType(const std::string &at1,
624	>	const std::string &at2,
625	>	const std::string &at3){
626	>	std::vector<std::string> keys;
627	>	keys.push_back(at1);
628	>	keys.push_back(at2);
629	>	keys.push_back(at3);
630	>	return bendTypeCont_.find(keys);
631	>	}
632	>
633	>	TorsionType* ForceField::getExactTorsionType(const std::string &at1,
634	>	const std::string &at2,
635	>	const std::string &at3,
636	>	const std::string &at4){
637	>	std::vector<std::string> keys;
638	>	keys.push_back(at1);
639	>	keys.push_back(at2);
640	>	keys.push_back(at3);
641	>	keys.push_back(at4);
642	>	return torsionTypeCont_.find(keys);
643	>	}
644	>
645	>	InversionType* ForceField::getExactInversionType(const std::string &at1,
646	>	const std::string &at2,
647	>	const std::string &at3,
648	>	const std::string &at4){
649	>	std::vector<std::string> keys;
650	>	keys.push_back(at1);
651	>	keys.push_back(at2);
652	>	keys.push_back(at3);
653	>	keys.push_back(at4);
654	>	return inversionTypeCont_.find(keys);
655	>	}
656	>
657	>	NonBondedInteractionType* ForceField::getExactNonBondedInteractionType(const std::string &at1, const std::string &at2){
658	>	std::vector<std::string> keys;
659	>	keys.push_back(at1);
660	>	keys.push_back(at2);
661	>	return nonBondedInteractionTypeCont_.find(keys);
662	>	}
663	>
664	>
665	>	bool ForceField::addAtomType(const std::string &at, AtomType* atomType) {
666	>	std::vector<std::string> keys;
667	>	keys.push_back(at);
668	>	atypeIdentToName[atomType->getIdent()] = at;
669	>	return atomTypeCont_.add(keys, atomType);
670	>	}
671	>
672	>	bool ForceField::replaceAtomType(const std::string &at, AtomType* atomType) {
673	>	std::vector<std::string> keys;
674	>	keys.push_back(at);
675	>	atypeIdentToName[atomType->getIdent()] = at;
676	>	return atomTypeCont_.replace(keys, atomType);
677	>	}
678	>
679	>	bool ForceField::addBondType(const std::string &at1, const std::string &at2,
680	>	BondType* bondType) {
681	>	std::vector<std::string> keys;
682	>	keys.push_back(at1);
683	>	keys.push_back(at2);
684	>	return bondTypeCont_.add(keys, bondType);
685	>	}
686	>
687	>	bool ForceField::addBendType(const std::string &at1, const std::string &at2,
688	>	const std::string &at3, BendType* bendType) {
689	>	std::vector<std::string> keys;
690	>	keys.push_back(at1);
691	>	keys.push_back(at2);
692	>	keys.push_back(at3);
693	>	return bendTypeCont_.add(keys, bendType);
694	>	}
695	>
696	>	bool ForceField::addTorsionType(const std::string &at1,
697	>	const std::string &at2,
698	>	const std::string &at3,
699	>	const std::string &at4,
700	>	TorsionType* torsionType) {
701	>	std::vector<std::string> keys;
702	>	keys.push_back(at1);
703	>	keys.push_back(at2);
704	>	keys.push_back(at3);
705	>	keys.push_back(at4);
706	>	return torsionTypeCont_.add(keys, torsionType);
707	>	}
708	>
709	>	bool ForceField::addInversionType(const std::string &at1,
710	>	const std::string &at2,
711	>	const std::string &at3,
712	>	const std::string &at4,
713	>	InversionType* inversionType) {
714	>	std::vector<std::string> keys;
715	>	keys.push_back(at1);
716	>	keys.push_back(at2);
717	>	keys.push_back(at3);
718	>	keys.push_back(at4);
719	>	return inversionTypeCont_.add(keys, inversionType);
720	>	}
721	>
722	>	bool ForceField::addNonBondedInteractionType(const std::string &at1,
723	>	const std::string &at2,
724	>	NonBondedInteractionType* nbiType) {
725	>	std::vector<std::string> keys;
726	>	keys.push_back(at1);
727	>	keys.push_back(at2);
728	>	return nonBondedInteractionTypeCont_.add(keys, nbiType);
729	>	}
730	>
731	>	RealType ForceField::getRcutFromAtomType(AtomType* at) {
732	>	RealType rcut(0.0);
733	>
734	>	LennardJonesAdapter lja = LennardJonesAdapter(at);
735	>	if (lja.isLennardJones()) {
736	>	rcut = 2.5 * lja.getSigma();
737	>	}
738	>	EAMAdapter ea = EAMAdapter(at);
739	>	if (ea.isEAM()) {
740	>	rcut = max(rcut, ea.getRcut());
741	>	}
742	>	SuttonChenAdapter sca = SuttonChenAdapter(at);
743	>	if (sca.isSuttonChen()) {
744	>	rcut = max(rcut, 2.0 * sca.getAlpha());
745	>	}
746	>	GayBerneAdapter gba = GayBerneAdapter(at);
747	>	if (gba.isGayBerne()) {
748	>	rcut = max(rcut, 2.5 * sqrt(2.0) * max(gba.getD(), gba.getL()));
749	>	}
750	>	StickyAdapter sa = StickyAdapter(at);
751	>	if (sa.isSticky()) {
752	>	rcut = max(rcut, max(sa.getRu(), sa.getRup()));
753	>	}
754	>
755	>	return rcut;
756	>	}
757	>
758	>
759	>	ifstrstream* ForceField::openForceFieldFile(const std::string& filename) {
760	>	std::string forceFieldFilename(filename);
761	>	ifstrstream* ffStream = new ifstrstream();
762	>
763	>	//try to open the force filed file in current directory first
764	>	ffStream->open(forceFieldFilename.c_str());
765	>	if(!ffStream->is_open()){
766	>
767	>	forceFieldFilename = ffPath_ + "/" + forceFieldFilename;
768	>	ffStream->open( forceFieldFilename.c_str() );
769	>
770	>	//if current directory does not contain the force field file,
771	>	//try to open it in the path
772	>	if(!ffStream->is_open()){
773	>
774	>	sprintf( painCave.errMsg,
775	>	"Error opening the force field parameter file:\n"
776	>	"\t%s\n"
777	>	"\tHave you tried setting the FORCE_PARAM_PATH environment "
778	>	"variable?\n",
779	>	forceFieldFilename.c_str() );
780	>	painCave.severity = OPENMD_ERROR;
781	>	painCave.isFatal = 1;
782	>	simError();
783	>	}
784	>	}
785	>	return ffStream;
786	>	}
787	>
788	>	} //end namespace OpenMD

Comparing:
trunk/src/UseTheForce/ForceField.cpp (property svn:keywords), Revision 206 by gezelter, Thu Nov 4 20:51:23 2004 UTC vs.
trunk/src/brains/ForceField.cpp (property svn:keywords), Revision 1880 by gezelter, Mon Jun 17 18:28:30 2013 UTC

#	Line 0 | Line 1
1	+	Author Id Revision Date

Diff Legend

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