1 // -*- C++ -*-
2 // KimParameterProvider.cpp: Get parameters from a string (for OpenKIM).
3 //
4 // Copyright (C) 2012-2013 Jakob Schiotz and the Department of Physics,
5 // Technical University of Denmark. Email: schiotz@fysik.dtu.dk
6 //
7 // This file is part of Asap version 3.
8 // Asap is released under the GNU Lesser Public License (LGPL) version 3.
9 // However, the parts of Asap distributed within the OpenKIM project
10 // (including this file) are also released under the Common Development
11 // and Distribution License (CDDL) version 1.0.
12 //
13 // This program is free software: you can redistribute it and/or
14 // modify it under the terms of the GNU Lesser General Public License
15 // version 3 as published by the Free Software Foundation. Permission
16 // to use other versions of the GNU Lesser General Public License may
17 // granted by Jakob Schiotz or the head of department of the
18 // Department of Physics, Technical University of Denmark, as
19 // described in section 14 of the GNU General Public License.
20 //
21 // This program is distributed in the hope that it will be useful,
22 // but WITHOUT ANY WARRANTY; without even the implied warranty of
23 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 // GNU General Public License for more details.
25 //
26 // You should have received a copy of the GNU General Public License
27 // and the GNU Lesser Public License along with this program. If not,
28 // see <http://www.gnu.org/licenses/>.
29
30
31 #include "KimParameterProvider.h"
32 #include "Asap.h"
33 #include "Exception.h"
34 // #define ASAPDEBUG
35 #include "Debug.h"
36 #include <iostream>
37 #include <fstream>
38 #include <cmath>
39 #include <set>
40 #include <stdlib.h>
41
KimParameterProvider(KIM::ModelDriverCreate * const modelDriverCreate,std::vector<std::string> & parameter_filenames,KIM::LengthUnit const requestedLengthUnit,KIM::EnergyUnit const requestedEnergyUnit,KIM::ChargeUnit const requestedChargeUnit,KIM::TemperatureUnit const requestedTemperatureUnit,KIM::TimeUnit const requestedTimeUnit)42 KimParameterProvider::KimParameterProvider(KIM::ModelDriverCreate *const modelDriverCreate,
43 std::vector<std::string> ¶meter_filenames,
44 KIM::LengthUnit const requestedLengthUnit,
45 KIM::EnergyUnit const requestedEnergyUnit,
46 KIM::ChargeUnit const requestedChargeUnit,
47 KIM::TemperatureUnit const requestedTemperatureUnit,
48 KIM::TimeUnit const requestedTimeUnit)
49 {
50 DEBUGPRINT;
51 int error;
52
53 // Register units
54 error = modelDriverCreate->SetUnits(requestedLengthUnit,
55 requestedEnergyUnit,
56 KIM::CHARGE_UNIT::unused,
57 KIM::TEMPERATURE_UNIT::unused,
58 KIM::TIME_UNIT::unused);
59 if (error)
60 throw AsapError("Unable to set units to requested values");
61
62 // define default base units
63 KIM::LengthUnit fromLength = KIM::LENGTH_UNIT::A;
64 KIM::EnergyUnit fromEnergy = KIM::ENERGY_UNIT::eV;
65 KIM::ChargeUnit fromCharge = KIM::CHARGE_UNIT::e;
66 KIM::TemperatureUnit fromTemperature = KIM::TEMPERATURE_UNIT::K;
67 KIM::TimeUnit fromTime = KIM::TIME_UNIT::ps;
68
69 // Prepare to convert all parameters to other units.
70 double energy_conv = 1.0;
71 error = modelDriverCreate->ConvertUnit(fromLength,
72 fromEnergy,
73 fromCharge,
74 fromTemperature,
75 fromTime,
76 requestedLengthUnit,
77 requestedEnergyUnit,
78 requestedChargeUnit,
79 requestedTemperatureUnit,
80 requestedTimeUnit,
81 0.0,
82 1.0,
83 0.0,
84 0.0,
85 0.0,
86 &energy_conv);
87 if (error)
88 throw AsapError("Unable to convert energy unit");
89
90 double len_conv = 1.0;
91 error = modelDriverCreate->ConvertUnit(fromLength,
92 fromEnergy,
93 fromCharge,
94 fromTemperature,
95 fromTime,
96 requestedLengthUnit,
97 requestedEnergyUnit,
98 requestedChargeUnit,
99 requestedTemperatureUnit,
100 requestedTimeUnit,
101 1.0,
102 0.0,
103 0.0,
104 0.0,
105 0.0,
106 &len_conv);
107 if (error)
108 throw AsapError("Unable to convert length unit");
109
110 double inv_len_conv = 1.0;
111 error = modelDriverCreate->ConvertUnit(fromLength,
112 fromEnergy,
113 fromCharge,
114 fromTemperature,
115 fromTime,
116 requestedLengthUnit,
117 requestedEnergyUnit,
118 requestedChargeUnit,
119 requestedTemperatureUnit,
120 requestedTimeUnit,
121 -1.0,
122 0.0,
123 0.0,
124 0.0,
125 0.0,
126 &inv_len_conv);
127 if (error)
128 throw AsapError("Unable to convert inverse length unit");
129
130 double inv_vol_conv = 1.0;
131 error = modelDriverCreate->ConvertUnit(fromLength,
132 fromEnergy,
133 fromCharge,
134 fromTemperature,
135 fromTime,
136 requestedLengthUnit,
137 requestedEnergyUnit,
138 requestedChargeUnit,
139 requestedTemperatureUnit,
140 requestedTimeUnit,
141 -3.0,
142 0.0,
143 0.0,
144 0.0,
145 0.0,
146 &inv_vol_conv);
147 if (error)
148 throw AsapError("Unable to convert inverse length unit");
149
150 // Read the parameter file.
151 if (parameter_filenames.size() != 1)
152 throw AsapError("Expected a single parameter file, got ") << parameter_filenames.size();
153 std::ifstream pstr(parameter_filenames[0].c_str());
154 int numread = -1;
155 emt_parameters *p = NULL;
156 while(true)
157 {
158 std::string word;
159 double value;
160 pstr >> word >> value;
161 //std::cerr << "read(" << word << ")(" << value << ")" << std::endl;
162 numread++;
163 if (word == "NEWELEMENT")
164 {
165 std::string name;
166 pstr >> name;
167 numread = 0;
168 p = new emt_parameters;
169 p->Z = (int) value;
170 p->gamma1 = 0.0; // These are calculated later!
171 p->gamma2 = 0.0;
172 p->name = name;
173 KIM::SpeciesName const specName(name);
174 error = modelDriverCreate->SetSpeciesCode(name, value);
175 if (error)
176 throw AsapError("Cannot set species code ") << name << " to " << value;
177 }
178 else if (word == "V0")
179 p->V0 = value * energy_conv;
180 else if (word == "kappa")
181 p->kappa = value * inv_len_conv;
182 else if (word == "eta2")
183 p->eta2 = value * inv_len_conv;
184 else if (word == "n0")
185 p->neq = value * inv_vol_conv;
186 else if (word == "S0")
187 p->seq = value * len_conv;
188 else if (word == "E0")
189 p->e0 = value * energy_conv;
190 else if (word == "lambda")
191 p->lambda = value * inv_len_conv;
192 else if (word == "mass")
193 {
194 p->mass = value; // Never really used.
195 p->invmass = 1.0/value;
196 }
197 else if (word == "ENDELEMENT")
198 {
199 if (numread != 9)
200 {
201 std::cerr << "Wrong number of parameters for element " << p->name << std::endl;
202 throw AsapError("Error in parameterfile");
203 }
204 p->index = params.size(); // Count the elements, starting at zero.
205 params.push_back(p);
206 p = NULL;
207 }
208 else if (word == "STOP")
209 break;
210 else
211 {
212 std::cerr << "Unknown keyword in parameter file: " << word << std::endl;
213 throw AsapError("Error in parameterfile");
214 }
215 }
216 }
217
GetNewParameters(int element)218 emt_parameters *KimParameterProvider::GetNewParameters(int element)
219 {
220 throw AsapError("Element not supported (missing from parameter file): number = ") << element;
221 }
222
223