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 475 by tim, Tue Apr 12 18:30:37 2005 UTC vs.
branches/development/src/brains/ForceField.cpp (file contents), Revision 1725 by gezelter, Sat May 26 18:13:43 2012 UTC

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

Comparing:
trunk/src/UseTheForce/ForceField.cpp (property svn:keywords), Revision 475 by tim, Tue Apr 12 18:30:37 2005 UTC vs.
branches/development/src/brains/ForceField.cpp (property svn:keywords), Revision 1725 by gezelter, Sat May 26 18:13:43 2012 UTC

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

Diff Legend

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