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 1277 by gezelter, Mon Jul 14 12:35:58 2008 UTC vs.
trunk/src/brains/ForceField.cpp (file contents), Revision 1790 by gezelter, Thu Aug 30 17:18:22 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		/**
#	Line 48 | Line 49
49		*/
50
51		#include <algorithm>
52	<	#include "UseTheForce/ForceField.hpp"
52	>	#include "brains/ForceField.hpp"
53		#include "utils/simError.h"
53	–	#include "utils/Tuple.hpp"
54	–	#include "UseTheForce/DarkSide/atype_interface.h"
55	–	#include "UseTheForce/DarkSide/fForceOptions_interface.h"
56	–	#include "UseTheForce/DarkSide/switcheroo_interface.h"
57	–	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 );
92		} else {
93		ffPath_ = tempPath;
94		}
95	+
96	+	setForceFieldFileName(ffName + ".frc");
97	+
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	+	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	<	ForceField::~ForceField() {
149	<	deleteAtypes();
150	<	deleteSwitch();
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	+	}
185	+
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,
#	Line 124 | Line 238 | namespace oopse {
238		myKeys.push_back((*i)->getName());
239		myKeys.push_back((*j)->getName());
240
127	–	std::cerr << "looking for " << myKeys[0] << " " << myKeys[1] << "\n";
241		BondType* bondType = bondTypeCont_.find(myKeys);
242		if (bondType) {
243		foundBonds.push_back(std::make_pair(bondTypeScore, myKeys));
#	Line 142 | Line 255 | namespace oopse {
255		int bestScore = foundBonds[0].first;
256		std::vector<std::string> theKeys = foundBonds[0].second;
257
145	–	std::cout << "best matching bond = " << theKeys[0] << "\t" << theKeys[1]  << "\t(score = "<< bestScore << ")\n";
258		BondType* bestType = bondTypeCont_.find(theKeys);
259
260		return bestType;
#	Line 226 | Line 338 | namespace oopse {
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;
341	>	std::vector<std::string> theKeys = foundBends[0].third;
342
231	–	std::cout << "best matching bend = " << theKeys[0] << "\t" <<theKeys[1]  << "\t" << theKeys[2] << "\t(scores = "<< jscore << "\t" << ikscore << ")\n";
232	–
343		BendType* bestType = bendTypeCont_.find(theKeys);
344		return bestType;
345		} else {
#	Line 332 | Line 442 | namespace oopse {
442		int ilscore = foundTorsions[0].second;
443		std::vector<std::string> theKeys = foundTorsions[0].third;
444
335	–	std::cout << "best matching torsion = " << theKeys[0] << "\t" <<theKeys[1]  << "\t" << theKeys[2] << "\t" << theKeys[3] << "\t(scores = "<< jkscore << "\t" << ilscore << ")\n";
336	–
445		TorsionType* bestType = torsionTypeCont_.find(theKeys);
446		return bestType;
447		} else {
#	Line 354 | Line 462 | namespace oopse {
462		keys.push_back(at4);
463
464		//try exact match first
465	<	InversionType* inversionType = inversionTypeCont_.find(keys);
465	>	InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(keys);
466		if (inversionType) {
467		return inversionType;
468		} else {
#	Line 416 | Line 524 | namespace oopse {
524		myKeys.push_back((*k)->getName());
525		myKeys.push_back((*l)->getName());
526
527	<	InversionType* inversionType = inversionTypeCont_.find(myKeys);
527	>	InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(myKeys);
528		if (inversionType) {
529		foundInversions.push_back( make_tuple3(Iscore, JKLscore, myKeys) );
530		}
#	Line 435 | Line 543 | namespace oopse {
543		int jklscore = foundInversions[0].second;
544		std::vector<std::string> theKeys = foundInversions[0].third;
545
546	<	std::cout << "best matching inversion = " << theKeys[0] << "\t" <<theKeys[1]  << "\t" << theKeys[2] << "\t" << theKeys[3] << "\t(scores = "<< iscore << "\t" << jklscore << ")\n";
439	<
440	<	InversionType* bestType = inversionTypeCont_.find(theKeys);
546	>	InversionType* bestType = inversionTypeCont_.permutedFindSkippingFirstElement(theKeys);
547		return bestType;
548		} else {
549		//if no exact match found, try wild card match
#	Line 447 | Line 553 | namespace oopse {
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);
#	Line 456 | Line 563 | namespace oopse {
563		if (nbiType) {
564		return nbiType;
565		} else {
566	<	//if no exact match found, try wild card match
567	<	return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_);
568	<	}
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,
#	Line 514 | Line 675 | namespace oopse {
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		}
681
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;
#	Line 570 | Line 739 | namespace oopse {
739		}
740
741		RealType ForceField::getRcutFromAtomType(AtomType* at) {
742	<	/**@todo */
574	<	GenericData* data;
575	<	RealType rcut = 0.0;
742	>	RealType rcut(0.0);
743
744	<	if (at->isLennardJones()) {
745	<	data = at->getPropertyByName("LennardJones");
746	<	if (data != NULL) {
580	<	LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
581	<
582	<	if (ljData != NULL) {
583	<	LJParam ljParam = ljData->getData();
584	<
585	<	//by default use 2.5*sigma as cutoff radius
586	<	rcut = 2.5 * ljParam.sigma;
587	<
588	<	} else {
589	<	sprintf( painCave.errMsg,
590	<	"Can not cast GenericData to LJParam\n");
591	<	painCave.severity = OOPSE_ERROR;
592	<	painCave.isFatal = 1;
593	<	simError();
594	<	}
595	<	} else {
596	<	sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n");
597	<	painCave.severity = OOPSE_ERROR;
598	<	painCave.isFatal = 1;
599	<	simError();
600	<	}
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		}
767
#	Line 608 | Line 771 | namespace oopse {
771		ifstrstream* ffStream = new ifstrstream();
772
773		//try to open the force filed file in current directory first
774	<	ffStream->open(forceFieldFilename.c_str());
774	>	ffStream->open(forceFieldFilename.c_str(), ifstream::in | ifstream::binary);
775	>
776		if(!ffStream->is_open()){
777
778		forceFieldFilename = ffPath_ + "/" + forceFieldFilename;
779	<	ffStream->open( forceFieldFilename.c_str() );
779	>	ffStream->open( forceFieldFilename.c_str(),
780	>	ifstream::in | ifstream::binary );
781
782		//if current directory does not contain the force field file,
783		//try to open it in the path
#	Line 624 | Line 789 | namespace oopse {
789		"\tHave you tried setting the FORCE_PARAM_PATH environment "
790		"variable?\n",
791		forceFieldFilename.c_str() );
792	<	painCave.severity = OOPSE_ERROR;
792	>	painCave.severity = OPENMD_ERROR;
793		painCave.isFatal = 1;
794		simError();
795		}
#	Line 632 | Line 797 | namespace oopse {
797		return ffStream;
798		}
799
800	<	void ForceField::setFortranForceOptions(){
636	<	ForceOptions theseFortranOptions;
637	<	forceFieldOptions_.makeFortranOptions(theseFortranOptions);
638	<	setfForceOptions(&theseFortranOptions);
639	<	}
640	<	} //end namespace oopse
800	>	} //end namespace OpenMD

Comparing:
trunk/src/UseTheForce/ForceField.cpp (property svn:keywords), Revision 1277 by gezelter, Mon Jul 14 12:35:58 2008 UTC vs.
trunk/src/brains/ForceField.cpp (property svn:keywords), Revision 1790 by gezelter, Thu Aug 30 17:18:22 2012 UTC

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

Diff Legend

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