ViewVC Help

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

root/OpenMD/branches/development/src/brains/ForceField.cpp

Comparing:
trunk/src/UseTheForce/ForceField.cpp (file contents), Revision 939 by gezelter, Thu Apr 20 18:24:24 2006 UTC vs.
branches/development/src/brains/ForceField.cpp (file contents), Revision 1808 by gezelter, Mon Oct 22 20:42:10 2012 UTC

#	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
46	–	* @time 22:51am
47		* @version 1.0
48		*/
49
50	<	#include "UseTheForce/ForceField.hpp"
50	>	#include <algorithm>
51	>	#include "brains/ForceField.hpp"
52		#include "utils/simError.h"
52	–	#include "UseTheForce/DarkSide/atype_interface.h"
53	–	#include "UseTheForce/DarkSide/fForceOptions_interface.h"
54	–	#include "UseTheForce/DarkSide/switcheroo_interface.h"
55	–	namespace oopse {
53
54	<	ForceField::ForceField() {
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	>	#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	>	namespace OpenMD {
82	>
83	>	ForceField::ForceField(std::string ffName) {
84	>
85		char* tempPath;
86		tempPath = getenv("FORCE_PARAM_PATH");
87	<
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	+	setForceFieldFileName(ffName + ".frc");
96	+
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	+	void ForceField::parse(const std::string& filename) {
145	+	ifstrstream* ffStream;
146
147	<	ForceField::~ForceField() {
148	<	deleteAtypes();
149	<	deleteSwitch();
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	+	/**
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	<	BondType* ForceField::getBondType(const std::string &at1, const std::string &at2) {
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);
#	Line 88 | Line 204 | namespace oopse {
204		if (bondType) {
205		return bondType;
206		} else {
207	<	//if no exact match found, try wild card match
208	<	return bondTypeCont_.find(keys, wildCardAtomTypeName_);
209	<	}
207	>	AtomType* atype1;
208	>	AtomType* atype2;
209	>	std::vector<std::string> at1key;
210	>	at1key.push_back(at1);
211	>	atype1 = atomTypeCont_.find(at1key);
212	>
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	<	BendType* ForceField::getBendType(const std::string &at1, const std::string &at2,
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.size() > 0) {
251	>	// sort the foundBonds by the score:
252	>	std::sort(foundBonds.begin(), foundBonds.end());
253	>
254	>	int bestScore = foundBonds[0].first;
255	>	std::vector<std::string> theKeys = foundBonds[0].second;
256	>
257	>	BondType* bestType = bondTypeCont_.find(theKeys);
258	>
259	>	return bestType;
260	>	} else {
261	>	//if no exact match found, try wild card match
262	>	return bondTypeCont_.find(keys, wildCardAtomTypeName_);
263	>	}
264	>	}
265	>	}
266	>
267	>	BendType* ForceField::getBendType(const std::string &at1,
268	>	const std::string &at2,
269		const std::string &at3) {
270		std::vector<std::string> keys;
271		keys.push_back(at1);
#	Line 106 | Line 277 | namespace oopse {
277		if (bendType) {
278		return bendType;
279		} else {
280	<	//if no exact match found, try wild card match
281	<	return bendTypeCont_.find(keys, wildCardAtomTypeName_);
280	>
281	>	AtomType* atype1;
282	>	AtomType* atype2;
283	>	AtomType* atype3;
284	>	std::vector<std::string> at1key;
285	>	at1key.push_back(at1);
286	>	atype1 = atomTypeCont_.find(at1key);
287	>
288	>	std::vector<std::string> at2key;
289	>	at2key.push_back(at2);
290	>	atype2 = atomTypeCont_.find(at2key);
291	>
292	>	std::vector<std::string> at3key;
293	>	at3key.push_back(at3);
294	>	atype3 = atomTypeCont_.find(at3key);
295	>
296	>	// query atom types for their chains of responsibility
297	>	std::vector<AtomType*> at1Chain = atype1->allYourBase();
298	>	std::vector<AtomType*> at2Chain = atype2->allYourBase();
299	>	std::vector<AtomType*> at3Chain = atype3->allYourBase();
300	>
301	>	std::vector<AtomType*>::iterator i;
302	>	std::vector<AtomType*>::iterator j;
303	>	std::vector<AtomType*>::iterator k;
304	>
305	>	int ii = 0;
306	>	int jj = 0;
307	>	int kk = 0;
308	>	int IKscore;
309	>
310	>	std::vector<tuple3<int, int, std::vector<std::string> > > foundBends;
311	>
312	>	for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
313	>	ii = 0;
314	>	for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
315	>	kk = 0;
316	>	for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
317	>
318	>	IKscore = ii + kk;
319	>
320	>	std::vector<std::string> myKeys;
321	>	myKeys.push_back((*i)->getName());
322	>	myKeys.push_back((*j)->getName());
323	>	myKeys.push_back((*k)->getName());
324	>
325	>	BendType* bendType = bendTypeCont_.find(myKeys);
326	>	if (bendType) {
327	>	foundBends.push_back( make_tuple3(jj, IKscore, myKeys) );
328	>	}
329	>	kk++;
330	>	}
331	>	ii++;
332	>	}
333	>	jj++;
334	>	}
335	>
336	>	if (foundBends.size() > 0) {
337	>	std::sort(foundBends.begin(), foundBends.end());
338	>	int jscore = foundBends[0].first;
339	>	int ikscore = foundBends[0].second;
340	>	std::vector<std::string> theKeys = foundBends[0].third;
341	>
342	>	BendType* bestType = bendTypeCont_.find(theKeys);
343	>	return bestType;
344	>	} else {
345	>	//if no exact match found, try wild card match
346	>	return bendTypeCont_.find(keys, wildCardAtomTypeName_);
347	>	}
348		}
349		}
350
351	<	TorsionType* ForceField::getTorsionType(const std::string &at1, const std::string &at2,
352	<	const std::string &at3, const std::string &at4) {
351	>	TorsionType* ForceField::getTorsionType(const std::string &at1,
352	>	const std::string &at2,
353	>	const std::string &at3,
354	>	const std::string &at4) {
355		std::vector<std::string> keys;
356		keys.push_back(at1);
357		keys.push_back(at2);
358		keys.push_back(at3);
359		keys.push_back(at4);
360
361	+
362	+	//try exact match first
363		TorsionType* torsionType = torsionTypeCont_.find(keys);
364		if (torsionType) {
365		return torsionType;
366		} else {
367	<	//if no exact match found, try wild card match
368	<	return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
369	<	}
370	<
371	<	return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
367	>
368	>	AtomType* atype1;
369	>	AtomType* atype2;
370	>	AtomType* atype3;
371	>	AtomType* atype4;
372	>	std::vector<std::string> at1key;
373	>	at1key.push_back(at1);
374	>	atype1 = atomTypeCont_.find(at1key);
375	>
376	>	std::vector<std::string> at2key;
377	>	at2key.push_back(at2);
378	>	atype2 = atomTypeCont_.find(at2key);
379	>
380	>	std::vector<std::string> at3key;
381	>	at3key.push_back(at3);
382	>	atype3 = atomTypeCont_.find(at3key);
383	>
384	>	std::vector<std::string> at4key;
385	>	at4key.push_back(at4);
386	>	atype4 = atomTypeCont_.find(at4key);
387	>
388	>	// query atom types for their chains of responsibility
389	>	std::vector<AtomType*> at1Chain = atype1->allYourBase();
390	>	std::vector<AtomType*> at2Chain = atype2->allYourBase();
391	>	std::vector<AtomType*> at3Chain = atype3->allYourBase();
392	>	std::vector<AtomType*> at4Chain = atype4->allYourBase();
393	>
394	>	std::vector<AtomType*>::iterator i;
395	>	std::vector<AtomType*>::iterator j;
396	>	std::vector<AtomType*>::iterator k;
397	>	std::vector<AtomType*>::iterator l;
398	>
399	>	int ii = 0;
400	>	int jj = 0;
401	>	int kk = 0;
402	>	int ll = 0;
403	>	int ILscore;
404	>	int JKscore;
405	>
406	>	std::vector<tuple3<int, int, std::vector<std::string> > > foundTorsions;
407	>
408	>	for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
409	>	kk = 0;
410	>	for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
411	>	ii = 0;
412	>	for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
413	>	ll = 0;
414	>	for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
415	>
416	>	ILscore = ii + ll;
417	>	JKscore = jj + kk;
418
419	+	std::vector<std::string> myKeys;
420	+	myKeys.push_back((*i)->getName());
421	+	myKeys.push_back((*j)->getName());
422	+	myKeys.push_back((*k)->getName());
423	+	myKeys.push_back((*l)->getName());
424	+
425	+	TorsionType* torsionType = torsionTypeCont_.find(myKeys);
426	+	if (torsionType) {
427	+	foundTorsions.push_back( make_tuple3(JKscore, ILscore, myKeys) );
428	+	}
429	+	ll++;
430	+	}
431	+	ii++;
432	+	}
433	+	kk++;
434	+	}
435	+	jj++;
436	+	}
437	+
438	+	if (foundTorsions.size() > 0) {
439	+	std::sort(foundTorsions.begin(), foundTorsions.end());
440	+	int jkscore = foundTorsions[0].first;
441	+	int ilscore = foundTorsions[0].second;
442	+	std::vector<std::string> theKeys = foundTorsions[0].third;
443	+
444	+	TorsionType* bestType = torsionTypeCont_.find(theKeys);
445	+	return bestType;
446	+	} else {
447	+	//if no exact match found, try wild card match
448	+	return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
449	+	}
450	+	}
451		}
452
453	<	BondType* ForceField::getExactBondType(const std::string &at1, const std::string &at2){
453	>	InversionType* ForceField::getInversionType(const std::string &at1,
454	>	const std::string &at2,
455	>	const std::string &at3,
456	>	const std::string &at4) {
457	>	std::vector<std::string> keys;
458	>	keys.push_back(at1);
459	>	keys.push_back(at2);
460	>	keys.push_back(at3);
461	>	keys.push_back(at4);
462	>
463	>	//try exact match first
464	>	InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(keys);
465	>	if (inversionType) {
466	>	return inversionType;
467	>	} else {
468	>
469	>	AtomType* atype1;
470	>	AtomType* atype2;
471	>	AtomType* atype3;
472	>	AtomType* atype4;
473	>	std::vector<std::string> at1key;
474	>	at1key.push_back(at1);
475	>	atype1 = atomTypeCont_.find(at1key);
476	>
477	>	std::vector<std::string> at2key;
478	>	at2key.push_back(at2);
479	>	atype2 = atomTypeCont_.find(at2key);
480	>
481	>	std::vector<std::string> at3key;
482	>	at3key.push_back(at3);
483	>	atype3 = atomTypeCont_.find(at3key);
484	>
485	>	std::vector<std::string> at4key;
486	>	at4key.push_back(at4);
487	>	atype4 = atomTypeCont_.find(at4key);
488	>
489	>	// query atom types for their chains of responsibility
490	>	std::vector<AtomType*> at1Chain = atype1->allYourBase();
491	>	std::vector<AtomType*> at2Chain = atype2->allYourBase();
492	>	std::vector<AtomType*> at3Chain = atype3->allYourBase();
493	>	std::vector<AtomType*> at4Chain = atype4->allYourBase();
494	>
495	>	std::vector<AtomType*>::iterator i;
496	>	std::vector<AtomType*>::iterator j;
497	>	std::vector<AtomType*>::iterator k;
498	>	std::vector<AtomType*>::iterator l;
499	>
500	>	int ii = 0;
501	>	int jj = 0;
502	>	int kk = 0;
503	>	int ll = 0;
504	>	int Iscore;
505	>	int JKLscore;
506	>
507	>	std::vector<tuple3<int, int, std::vector<std::string> > > foundInversions;
508	>
509	>	for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
510	>	kk = 0;
511	>	for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
512	>	ii = 0;
513	>	for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
514	>	ll = 0;
515	>	for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
516	>
517	>	Iscore = ii;
518	>	JKLscore = jj + kk + ll;
519	>
520	>	std::vector<std::string> myKeys;
521	>	myKeys.push_back((*i)->getName());
522	>	myKeys.push_back((*j)->getName());
523	>	myKeys.push_back((*k)->getName());
524	>	myKeys.push_back((*l)->getName());
525	>
526	>	InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(myKeys);
527	>	if (inversionType) {
528	>	foundInversions.push_back( make_tuple3(Iscore, JKLscore, myKeys) );
529	>	}
530	>	ll++;
531	>	}
532	>	ii++;
533	>	}
534	>	kk++;
535	>	}
536	>	jj++;
537	>	}
538	>
539	>	if (foundInversions.size() > 0) {
540	>	std::sort(foundInversions.begin(), foundInversions.end());
541	>	int iscore = foundInversions[0].first;
542	>	int jklscore = foundInversions[0].second;
543	>	std::vector<std::string> theKeys = foundInversions[0].third;
544	>
545	>	InversionType* bestType = inversionTypeCont_.permutedFindSkippingFirstElement(theKeys);
546	>	return bestType;
547	>	} else {
548	>	//if no exact match found, try wild card match
549	>	return inversionTypeCont_.find(keys, wildCardAtomTypeName_);
550	>	}
551	>	}
552	>	}
553	>
554	>	NonBondedInteractionType* ForceField::getNonBondedInteractionType(const std::string &at1, const std::string &at2) {
555	>
556	>	std::vector<std::string> keys;
557	>	keys.push_back(at1);
558	>	keys.push_back(at2);
559	>
560	>	//try exact match first
561	>	NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(keys);
562	>	if (nbiType) {
563	>	return nbiType;
564	>	} else {
565	>	AtomType* atype1;
566	>	AtomType* atype2;
567	>	std::vector<std::string> at1key;
568	>	at1key.push_back(at1);
569	>	atype1 = atomTypeCont_.find(at1key);
570	>
571	>	std::vector<std::string> at2key;
572	>	at2key.push_back(at2);
573	>	atype2 = atomTypeCont_.find(at2key);
574	>
575	>	// query atom types for their chains of responsibility
576	>	std::vector<AtomType*> at1Chain = atype1->allYourBase();
577	>	std::vector<AtomType*> at2Chain = atype2->allYourBase();
578	>
579	>	std::vector<AtomType*>::iterator i;
580	>	std::vector<AtomType*>::iterator j;
581	>
582	>	int ii = 0;
583	>	int jj = 0;
584	>	int nbiTypeScore;
585	>
586	>	std::vector<std::pair<int, std::vector<std::string> > > foundNBI;
587	>
588	>	for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
589	>	jj = 0;
590	>	for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
591	>
592	>	nbiTypeScore = ii + jj;
593	>
594	>	std::vector<std::string> myKeys;
595	>	myKeys.push_back((*i)->getName());
596	>	myKeys.push_back((*j)->getName());
597	>
598	>	NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(myKeys);
599	>	if (nbiType) {
600	>	foundNBI.push_back(std::make_pair(nbiTypeScore, myKeys));
601	>	}
602	>	jj++;
603	>	}
604	>	ii++;
605	>	}
606	>
607	>
608	>	if (foundNBI.size() > 0) {
609	>	// sort the foundNBI by the score:
610	>	std::sort(foundNBI.begin(), foundNBI.end());
611	>
612	>	int bestScore = foundNBI[0].first;
613	>	std::vector<std::string> theKeys = foundNBI[0].second;
614	>
615	>	NonBondedInteractionType* bestType = nonBondedInteractionTypeCont_.find(theKeys);
616	>	return bestType;
617	>	} else {
618	>	//if no exact match found, try wild card match
619	>	return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_);
620	>	}
621	>	}
622	>	}
623	>
624	>	BondType* ForceField::getExactBondType(const std::string &at1,
625	>	const std::string &at2){
626		std::vector<std::string> keys;
627		keys.push_back(at1);
628		keys.push_back(at2);
629		return bondTypeCont_.find(keys);
630		}
631	<
632	<	BendType* ForceField::getExactBendType(const std::string &at1, const std::string &at2,
631	>
632	>	BendType* ForceField::getExactBendType(const std::string &at1,
633	>	const std::string &at2,
634		const std::string &at3){
635		std::vector<std::string> keys;
636		keys.push_back(at1);
#	Line 146 | Line 638 | namespace oopse {
638		keys.push_back(at3);
639		return bendTypeCont_.find(keys);
640		}
641	<
642	<	TorsionType* ForceField::getExactTorsionType(const std::string &at1, const std::string &at2,
643	<	const std::string &at3, const std::string &at4){
641	>
642	>	TorsionType* ForceField::getExactTorsionType(const std::string &at1,
643	>	const std::string &at2,
644	>	const std::string &at3,
645	>	const std::string &at4){
646		std::vector<std::string> keys;
647		keys.push_back(at1);
648		keys.push_back(at2);
#	Line 156 | Line 650 | namespace oopse {
650		keys.push_back(at4);
651		return torsionTypeCont_.find(keys);
652		}
653	+
654	+	InversionType* ForceField::getExactInversionType(const std::string &at1,
655	+	const std::string &at2,
656	+	const std::string &at3,
657	+	const std::string &at4){
658	+	std::vector<std::string> keys;
659	+	keys.push_back(at1);
660	+	keys.push_back(at2);
661	+	keys.push_back(at3);
662	+	keys.push_back(at4);
663	+	return inversionTypeCont_.find(keys);
664	+	}
665	+
666	+	NonBondedInteractionType* ForceField::getExactNonBondedInteractionType(const std::string &at1, const std::string &at2){
667	+	std::vector<std::string> keys;
668	+	keys.push_back(at1);
669	+	keys.push_back(at2);
670	+	return nonBondedInteractionTypeCont_.find(keys);
671	+	}
672	+
673	+
674		bool ForceField::addAtomType(const std::string &at, AtomType* atomType) {
675		std::vector<std::string> keys;
676		keys.push_back(at);
677	+	atypeIdentToName[atomType->getIdent()] = at;
678		return atomTypeCont_.add(keys, atomType);
679		}
680
681	<	bool ForceField::addBondType(const std::string &at1, const std::string &at2, BondType* bondType) {
681	>	bool ForceField::replaceAtomType(const std::string &at, AtomType* atomType) {
682		std::vector<std::string> keys;
683	+	keys.push_back(at);
684	+	atypeIdentToName[atomType->getIdent()] = at;
685	+	return atomTypeCont_.replace(keys, atomType);
686	+	}
687	+
688	+	bool ForceField::addBondType(const std::string &at1, const std::string &at2,
689	+	BondType* bondType) {
690	+	std::vector<std::string> keys;
691		keys.push_back(at1);
692		keys.push_back(at2);
693	<	return bondTypeCont_.add(keys, bondType);
170	<
693	>	return bondTypeCont_.add(keys, bondType);
694		}
695	<
695	>
696		bool ForceField::addBendType(const std::string &at1, const std::string &at2,
697		const std::string &at3, BendType* bendType) {
698		std::vector<std::string> keys;
#	Line 178 | Line 701 | namespace oopse {
701		keys.push_back(at3);
702		return bendTypeCont_.add(keys, bendType);
703		}
704	<
705	<	bool ForceField::addTorsionType(const std::string &at1, const std::string &at2,
706	<	const std::string &at3, const std::string &at4, TorsionType* torsionType) {
704	>
705	>	bool ForceField::addTorsionType(const std::string &at1,
706	>	const std::string &at2,
707	>	const std::string &at3,
708	>	const std::string &at4,
709	>	TorsionType* torsionType) {
710		std::vector<std::string> keys;
711		keys.push_back(at1);
712		keys.push_back(at2);
#	Line 189 | Line 715 | namespace oopse {
715		return torsionTypeCont_.add(keys, torsionType);
716		}
717
718	<	double ForceField::getRcutFromAtomType(AtomType* at) {
719	<	/**@todo */
720	<	GenericData* data;
721	<	double rcut = 0.0;
722	<
723	<	if (at->isLennardJones()) {
724	<	data = at->getPropertyByName("LennardJones");
725	<	if (data != NULL) {
726	<	LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
727	<
728	<	if (ljData != NULL) {
729	<	LJParam ljParam = ljData->getData();
730	<
731	<	//by default use 2.5*sigma as cutoff radius
732	<	rcut = 2.5 * ljParam.sigma;
733	<
734	<	} else {
735	<	sprintf( painCave.errMsg,
736	<	"Can not cast GenericData to LJParam\n");
737	<	painCave.severity = OOPSE_ERROR;
738	<	painCave.isFatal = 1;
739	<	simError();
740	<	}
741	<	} else {
742	<	sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n");
743	<	painCave.severity = OOPSE_ERROR;
744	<	painCave.isFatal = 1;
745	<	simError();
220	<	}
718	>	bool ForceField::addInversionType(const std::string &at1,
719	>	const std::string &at2,
720	>	const std::string &at3,
721	>	const std::string &at4,
722	>	InversionType* inversionType) {
723	>	std::vector<std::string> keys;
724	>	keys.push_back(at1);
725	>	keys.push_back(at2);
726	>	keys.push_back(at3);
727	>	keys.push_back(at4);
728	>	return inversionTypeCont_.add(keys, inversionType);
729	>	}
730	>
731	>	bool ForceField::addNonBondedInteractionType(const std::string &at1,
732	>	const std::string &at2,
733	>	NonBondedInteractionType* nbiType) {
734	>	std::vector<std::string> keys;
735	>	keys.push_back(at1);
736	>	keys.push_back(at2);
737	>	return nonBondedInteractionTypeCont_.add(keys, nbiType);
738	>	}
739	>
740	>	RealType ForceField::getRcutFromAtomType(AtomType* at) {
741	>	RealType rcut(0.0);
742	>
743	>	LennardJonesAdapter lja = LennardJonesAdapter(at);
744	>	if (lja.isLennardJones()) {
745	>	rcut = 2.5 * lja.getSigma();
746		}
747	+	EAMAdapter ea = EAMAdapter(at);
748	+	if (ea.isEAM()) {
749	+	rcut = max(rcut, ea.getRcut());
750	+	}
751	+	SuttonChenAdapter sca = SuttonChenAdapter(at);
752	+	if (sca.isSuttonChen()) {
753	+	rcut = max(rcut, 2.0 * sca.getAlpha());
754	+	}
755	+	GayBerneAdapter gba = GayBerneAdapter(at);
756	+	if (gba.isGayBerne()) {
757	+	rcut = max(rcut, 2.5 * sqrt(2.0) * max(gba.getD(), gba.getL()));
758	+	}
759	+	StickyAdapter sa = StickyAdapter(at);
760	+	if (sa.isSticky()) {
761	+	rcut = max(rcut, max(sa.getRu(), sa.getRup()));
762	+	}
763
764		return rcut;
765		}
766	+
767
226	–
768		ifstrstream* ForceField::openForceFieldFile(const std::string& filename) {
769		std::string forceFieldFilename(filename);
770		ifstrstream* ffStream = new ifstrstream();
#	Line 245 | Line 786 | namespace oopse {
786		"\tHave you tried setting the FORCE_PARAM_PATH environment "
787		"variable?\n",
788		forceFieldFilename.c_str() );
789	<	painCave.severity = OOPSE_ERROR;
789	>	painCave.severity = OPENMD_ERROR;
790		painCave.isFatal = 1;
791		simError();
792		}
793		}
253	–
794		return ffStream;
255	–
795		}
796
797	<	void ForceField::setFortranForceOptions(){
259	<	ForceOptions theseFortranOptions;
260	<	forceFieldOptions_.makeFortranOptions(theseFortranOptions);
261	<	setfForceOptions(&theseFortranOptions);
262	<	}
263	<	} //end namespace oopse
797	>	} //end namespace OpenMD

Comparing:
trunk/src/UseTheForce/ForceField.cpp (property svn:keywords), Revision 939 by gezelter, Thu Apr 20 18:24:24 2006 UTC vs.
branches/development/src/brains/ForceField.cpp (property svn:keywords), Revision 1808 by gezelter, Mon Oct 22 20:42:10 2012 UTC

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

Diff Legend

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