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> &parameter_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