File: quadraticImproperTorsion.C

package info (click to toggle)
ball 1.5.0%2Bgit20180813.37fc53c-11
  • links: PTS, VCS
  • area: main
  • in suites: bookworm
  • size: 239,924 kB
  • sloc: cpp: 326,149; ansic: 4,208; python: 2,303; yacc: 1,778; lex: 1,099; xml: 958; sh: 322; javascript: 164; makefile: 88
file content (269 lines) | stat: -rw-r--r-- 7,084 bytes parent folder | download | duplicates (8)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
// -*- Mode: C++; tab-width: 2; -*-
// vi: set ts=2:
//

#include <BALL/MOLMEC/PARAMETER/quadraticImproperTorsion.h>
#include <BALL/MOLMEC/PARAMETER/forceFieldParameters.h>

using namespace std;

namespace BALL 
{

	QuadraticImproperTorsion::QuadraticImproperTorsion() 
		:	ParameterSection(),
			number_of_atom_types_(0),
			torsions_(),
			torsion_hash_map_()
	{
	}


	QuadraticImproperTorsion::QuadraticImproperTorsion
		(const QuadraticImproperTorsion& torsion) 
		:	ParameterSection(torsion),
			number_of_atom_types_(torsion.number_of_atom_types_),
			torsions_(torsion.torsions_),
			torsion_hash_map_(torsion.torsion_hash_map_)
	{
	}


	QuadraticImproperTorsion::~QuadraticImproperTorsion() 
	{
		clear();

		valid_ = false;
	}


	void QuadraticImproperTorsion::clear() 
	{
		number_of_atom_types_ = 0;
		torsions_.clear();
		torsion_hash_map_.clear();

		ParameterSection::clear();
	}


	const QuadraticImproperTorsion& QuadraticImproperTorsion::operator =
		(const QuadraticImproperTorsion& torsion) 
	{
		ParameterSection::operator = (torsion);

		number_of_atom_types_ = torsion.number_of_atom_types_;
		torsions_ = torsion.torsions_;
		torsion_hash_map_ = torsion.torsion_hash_map_;

		return *this;
	}


	bool QuadraticImproperTorsion::operator ==
		(const QuadraticImproperTorsion& torsion) const 
	{
		return (ParameterSection::operator == (torsion)
			&& (number_of_atom_types_ == torsion.number_of_atom_types_)
			&& (torsions_ == torsion.torsions_)
			&& (torsion_hash_map_ == torsion.torsion_hash_map_));
	}


	bool QuadraticImproperTorsion::extractSection
		(Parameters& parameters, const String& section_name) 
	{
		return ParameterSection::extractSection(parameters, section_name);
	}

	bool QuadraticImproperTorsion::extractSection
		(ForceFieldParameters& parameters, const String& section_name) 
	{
		// first clear the fields

		clear();

		// check whether the parameters are valid
		if (!parameters.isValid())
		{
			return false;
		}
		
		// extract the section information
		if (!ParameterSection::extractSection(parameters, section_name))
		{
			Log.error() << "Could not find section " << section_name << " in parameter file!" << endl;
			return false;
		}
		
		// check whether all variables we need are defined, terminate otherwise
		if (!hasVariable("phase") || !hasVariable("k"))
		{
			Log.error() << "QuadraticImproperTorsion section (" << section_name << ") needs columns phase and k!" << endl;
			return false;
		}

		// build a two dimensional array of the atom types
		// loop variable
		const AtomTypes&	atom_types = parameters.getAtomTypes();
		number_of_atom_types_ = atom_types.getNumberOfTypes();

		// clear all old torsions
		//	- torsions_ is a vector containing Values objects
		//  - torsion_hash_map_ hashes the product of the atom types (I + number_of_atom_types_ * J +...)
		//		to a the index in torsions_
		torsions_.clear();
		torsion_hash_map_.clear();
		
		// determine the units of the phase and potential wall
		// (if given in options)
		float factor_phase = 1.0;
		if (options.has("unit_phase"))
		{
			if (options["unit_phase"] == "deg")
			{	
				factor_phase = Constants::PI / 180.0;
			}
		}

		float factor_k = 1.0;
		if (options.has("unit_k"))
		{
			if (options["unit_k"] == "kcal/mol")
			{
				factor_k = Constants::JOULE_PER_CAL;
			}
		}

		Atom::Type		type_I;
		Atom::Type		type_J;
		Atom::Type		type_K;
		Atom::Type		type_L;
		String				key;
		String				fields[5];

		StringHashMap<Index>::Iterator it;
		for (it = section_entries_.begin(); it != section_entries_.end(); ++it)
		{
			key = it->first;
			if (key.split(fields, 4) == 4)
			{
				// determine all atom types
				type_I = atom_types.getType(fields[0]);
				type_J = atom_types.getType(fields[1]);
				type_K = atom_types.getType(fields[2]);
				type_L = atom_types.getType(fields[3]);

				// create a new torsion and store 
				// it in the vector of torsions
				Size array_idx = (Size)torsions_.size();	
				Values v;
				v.phase = getValue(key, "phase").toFloat() * factor_phase;
				v.k = getValue(key, "k").toFloat() * factor_k;
				torsions_.push_back(v);
						
				// insert the array index and the atom type key into the
				// hash map

				// calculate a unique number for each possible combination of
				// atom types and use it to hash the torsion parameters	
				Size index = type_I + type_J * number_of_atom_types_ 
									 + type_K * number_of_atom_types_ * number_of_atom_types_ 
									 + type_L * number_of_atom_types_ * number_of_atom_types_ * number_of_atom_types_;
				torsion_hash_map_.insert(pair<Size, Size>(index, array_idx));
			} 
			else 
			{
				Log.error() << "Could not interpret key " << key << endl;
			}
		}

		return true;
	}


	bool QuadraticImproperTorsion::hasParameters
		(Atom::Type I, Atom::Type J, Atom::Type K, Atom::Type L) const 
	{
		if ((I < 0) || ((Size)I >= number_of_atom_types_))
		{
			return false;
		}

		if ((J < 0) || ((Size)J >= number_of_atom_types_))
		{
			return false;
		}

		if ((K < 0) || ((Size)K >= number_of_atom_types_))
		{
			return false;
		}

		if ((L < 0) || ((Size)L >= number_of_atom_types_))
		{
			return false;
		}

		// calculate the key for this combination of atom types
		Size index = I + number_of_atom_types_ * J 
							 + K * number_of_atom_types_ * number_of_atom_types_
							 + L * number_of_atom_types_ * number_of_atom_types_ * number_of_atom_types_;

		// and look it up in the hash table
		bool result = torsion_hash_map_.has(index);

		// now, check for wildcards at the inner positions
		if (!result)
		{
			index = I + number_of_atom_types_ * Atom::ANY_TYPE 
							 + Atom::ANY_TYPE * number_of_atom_types_ * number_of_atom_types_
							 + L * number_of_atom_types_ * number_of_atom_types_ * number_of_atom_types_;
			result = torsion_hash_map_.has(index);
		}

		return result;
	}


	QuadraticImproperTorsion::Values QuadraticImproperTorsion::getParameters
		(Atom::Type I, Atom::Type J, Atom::Type K, Atom::Type L) const 
	{
		QuadraticImproperTorsion::Values parameters;
		assignParameters(parameters, I, J, K, L);
		return parameters;
	}


	bool QuadraticImproperTorsion::assignParameters
		(QuadraticImproperTorsion::Values& parameters,
		 Atom::Type I, Atom::Type J, Atom::Type K, Atom::Type L) const 
	{
		// calculate the key for this combination of atom types
		Size index = I + number_of_atom_types_ * J 
							 + K * number_of_atom_types_ * number_of_atom_types_
							 + L * number_of_atom_types_ * number_of_atom_types_ * number_of_atom_types_;

		// and look it up in the hash table
		bool result = torsion_hash_map_.has(index);

		// now, check for wildcards at the inner positions
		if (!result)
		{
			index = I + number_of_atom_types_ * Atom::ANY_TYPE 
							 + Atom::ANY_TYPE * number_of_atom_types_ * number_of_atom_types_
							 + L * number_of_atom_types_ * number_of_atom_types_ * number_of_atom_types_;
			result = torsion_hash_map_.has(index);
		}

		// and look it up in the hash table
		if (result)
		{
			parameters = torsions_[torsion_hash_map_[index]];
		}

		return result;
	}
	 

} // namespace BALL