n2p2 - A neural network potential package
ElementCabana_impl.h
Go to the documentation of this file.
1// n2p2 - A neural network potential package
2// Copyright (C) 2018 Andreas Singraber (University of Vienna)
3// Copyright (C) 2020 Saaketh Desai and Sam Reeve
4//
5// This program is free software: you can redistribute it and/or modify
6// it under the terms of the GNU General Public License as published by
7// the Free Software Foundation, either version 3 of the License, or
8// (at your option) any later version.
9//
10// This program is distributed in the hope that it will be useful,
11// but WITHOUT ANY WARRANTY; without even the implied warranty of
12// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13// GNU General Public License for more details.
14//
15// You should have received a copy of the GNU General Public License
16// along with this program. If not, see <https://www.gnu.org/licenses/>.
17
18#include "utility.h"
19
20#include <algorithm> // std::sort, std::min, std::max
21#include <cstdlib> // atoi
22#include <iostream> // std::cerr
23#include <limits> // std::numeric_limits
24#include <stdexcept> // std::runtime_error
25
26using namespace std;
27
28namespace nnp
29{
30
31ElementCabana::ElementCabana( size_t const _index )
32 : Element()
33{
34 index = _index;
35 atomicEnergyOffset = 0.0;
36 neuralNetwork = NULL;
37}
38
39ElementCabana::~ElementCabana() {}
40
41template <class t_SF, class h_t_int>
42void ElementCabana::addSymmetryFunction( string const &parameters,
43 vector<string> elementStrings, int attype,
44 t_SF SF, double convLength,
45 h_t_int h_numSFperElem )
46{
47 vector<string> args = split( reduce( parameters ) );
48 size_t type = (size_t)atoi( args.at( 1 ).c_str() );
49 const char *estring;
50 int el = 0;
51
52 vector<string> splitLine = split( reduce( parameters ) );
53 if ( type == 2 )
54 {
55 estring = splitLine.at( 0 ).c_str();
56 for ( size_t i = 0; i < elementStrings.size(); ++i )
57 {
58 if ( strcmp( elementStrings[i].c_str(), estring ) == 0 )
59 el = i;
60 }
61 SF( attype, h_numSFperElem( attype ), 0 ) = el; // ec
62 SF( attype, h_numSFperElem( attype ), 1 ) = type; // type
63
64 estring = splitLine.at( 2 ).c_str();
65 for ( size_t i = 0; i < elementStrings.size(); ++i )
66 {
67 if ( strcmp( elementStrings[i].c_str(), estring ) == 0 )
68 el = i;
69 }
70 SF( attype, h_numSFperElem( attype ), 2 ) = el; // e1
71 // set e2 to arbit high number for ease in creating groups
72 SF( attype, h_numSFperElem( attype ), 3 ) = 100000; // e2
73
74 SF( attype, h_numSFperElem( attype ), 4 ) =
75 atof( splitLine.at( 3 ).c_str() ) /
76 ( convLength * convLength ); // eta
77 SF( attype, h_numSFperElem( attype ), 8 ) =
78 atof( splitLine.at( 4 ).c_str() ) * convLength; // rs
79 SF( attype, h_numSFperElem( attype ), 7 ) =
80 atof( splitLine.at( 5 ).c_str() ) * convLength; // rc
81
82 SF( attype, h_numSFperElem( attype ), 13 ) = h_numSFperElem( attype );
83 h_numSFperElem( attype )++;
84 }
85
86 else if ( type == 3 )
87 {
88 if ( type != (size_t)atoi( splitLine.at( 1 ).c_str() ) )
89 throw runtime_error( "ERROR: Incorrect symmetry function type.\n" );
90 estring = splitLine.at( 0 ).c_str();
91 for ( size_t i = 0; i < elementStrings.size(); ++i )
92 {
93 if ( strcmp( elementStrings[i].c_str(), estring ) == 0 )
94 el = i;
95 }
96 SF( attype, h_numSFperElem( attype ), 0 ) = el; // ec
97 SF( attype, h_numSFperElem( attype ), 1 ) = type; // type
98
99 estring = splitLine.at( 2 ).c_str();
100 for ( size_t i = 0; i < elementStrings.size(); ++i )
101 {
102 if ( strcmp( elementStrings[i].c_str(), estring ) == 0 )
103 el = i;
104 }
105 SF( attype, h_numSFperElem( attype ), 2 ) = el; // e1
106
107 estring = splitLine.at( 3 ).c_str();
108 for ( size_t i = 0; i < elementStrings.size(); ++i )
109 {
110 if ( strcmp( elementStrings[i].c_str(), estring ) == 0 )
111 el = i;
112 }
113
114 SF( attype, h_numSFperElem( attype ), 3 ) = el; // e2
115 SF( attype, h_numSFperElem( attype ), 4 ) =
116 atof( splitLine.at( 4 ).c_str() ) /
117 ( convLength * convLength ); // eta
118 SF( attype, h_numSFperElem( attype ), 5 ) =
119 atof( splitLine.at( 5 ).c_str() ); // lambda
120 SF( attype, h_numSFperElem( attype ), 6 ) =
121 atof( splitLine.at( 6 ).c_str() ); // zeta
122 SF( attype, h_numSFperElem( attype ), 7 ) =
123 atof( splitLine.at( 7 ).c_str() ) * convLength; // rc
124 // Shift parameter is optional.
125 if ( splitLine.size() > 8 )
126 SF( attype, h_numSFperElem( attype ), 8 ) =
127 atof( splitLine.at( 8 ).c_str() ) * convLength; // rs
128
129 T_INT e1 = SF( attype, h_numSFperElem( attype ), 2 );
130 T_INT e2 = SF( attype, h_numSFperElem( attype ), 3 );
131 if ( e1 > e2 )
132 {
133 size_t tmp = e1;
134 e1 = e2;
135 e2 = tmp;
136 }
137 SF( attype, h_numSFperElem( attype ), 2 ) = e1;
138 SF( attype, h_numSFperElem( attype ), 3 ) = e2;
139
140 T_FLOAT zeta = SF( attype, h_numSFperElem( attype ), 6 );
141 T_INT zetaInt = round( zeta );
142 if ( fabs( zeta - zetaInt ) <= numeric_limits<double>::min() )
143 SF( attype, h_numSFperElem( attype ), 9 ) = 1;
144 else
145 SF( attype, h_numSFperElem( attype ), 9 ) = 0;
146
147 SF( attype, h_numSFperElem( attype ), 13 ) = h_numSFperElem( attype );
148 h_numSFperElem( attype )++;
149 }
150 // TODO: Add this later
151 else if ( type == 9 )
152 {
153 }
154 else if ( type == 12 )
155 {
156 }
157 else if ( type == 13 )
158 {
159 }
160 else
161 {
162 throw runtime_error( "ERROR: Unknown symmetry function type.\n" );
163 }
164
165 return;
166}
167
168template <class t_SF, class h_t_int>
169void ElementCabana::sortSymmetryFunctions( t_SF SF, h_t_int h_numSFperElem,
170 int attype )
171{
172 int size = h_numSFperElem( attype );
173 h_t_int h_SFsort( "SortSort", size );
174 for ( int i = 0; i < size; ++i )
175 h_SFsort( i ) = i;
176
177 // naive insertion sort
178 int i, j, tmp;
179 for ( i = 1; i < size; ++i )
180 {
181 j = i;
182 // explicit condition for sort
183 while ( j > 0 &&
184 compareSF( SF, attype, h_SFsort( j - 1 ), h_SFsort( j ) ) )
185 {
186 tmp = h_SFsort( j );
187 h_SFsort( j ) = h_SFsort( j - 1 );
188 h_SFsort( j - 1 ) = tmp;
189 --j;
190 }
191 }
192
193 int tmpindex;
194 for ( int i = 0; i < size; ++i )
195 {
196 SF( attype, i, 13 ) = h_SFsort( i );
197 tmpindex = SF( attype, i, 13 );
198 SF( attype, tmpindex, 14 ) = i;
199 }
200
201 return;
202}
203
204template <class t_SF>
205bool ElementCabana::compareSF( t_SF SF, int attype, int index1, int index2 )
206{
207 if ( SF( attype, index2, 0 ) < SF( attype, index1, 0 ) )
208 return true; // ec
209 else if ( SF( attype, index2, 0 ) > SF( attype, index1, 0 ) )
210 return false;
211
212 if ( SF( attype, index2, 1 ) < SF( attype, index1, 1 ) )
213 return true; // type
214 else if ( SF( attype, index2, 1 ) > SF( attype, index1, 1 ) )
215 return false;
216
217 if ( SF( attype, index2, 11 ) < SF( attype, index1, 11 ) )
218 return true; // cutofftype
219 else if ( SF( attype, index2, 11 ) > SF( attype, index1, 11 ) )
220 return false;
221
222 if ( SF( attype, index2, 12 ) < SF( attype, index1, 12 ) )
223 return true; // cutoffalpha
224 else if ( SF( attype, index2, 12 ) > SF( attype, index1, 12 ) )
225 return false;
226
227 if ( SF( attype, index2, 7 ) < SF( attype, index1, 7 ) )
228 return true; // rc
229 else if ( SF( attype, index2, 7 ) > SF( attype, index1, 7 ) )
230 return false;
231
232 if ( SF( attype, index2, 4 ) < SF( attype, index1, 4 ) )
233 return true; // eta
234 else if ( SF( attype, index2, 4 ) > SF( attype, index1, 4 ) )
235 return false;
236
237 if ( SF( attype, index2, 8 ) < SF( attype, index1, 8 ) )
238 return true; // rs
239 else if ( SF( attype, index2, 8 ) > SF( attype, index1, 8 ) )
240 return false;
241
242 if ( SF( attype, index2, 6 ) < SF( attype, index1, 6 ) )
243 return true; // zeta
244 else if ( SF( attype, index2, 6 ) > SF( attype, index1, 6 ) )
245 return false;
246
247 if ( SF( attype, index2, 5 ) < SF( attype, index1, 5 ) )
248 return true; // lambda
249 else if ( SF( attype, index2, 5 ) > SF( attype, index1, 5 ) )
250 return false;
251
252 if ( SF( attype, index2, 2 ) < SF( attype, index1, 2 ) )
253 return true; // e1
254 else if ( SF( attype, index2, 2 ) > SF( attype, index1, 2 ) )
255 return false;
256
257 if ( SF( attype, index2, 3 ) < SF( attype, index1, 3 ) )
258 return true; // e2
259 else if ( SF( attype, index2, 3 ) > SF( attype, index1, 3 ) )
260 return false;
261
262 else
263 return false;
264}
265
266template <class t_SF, class h_t_int>
267vector<string>
268ElementCabana::infoSymmetryFunctionParameters( t_SF SF, int attype,
269 h_t_int h_numSFperElem ) const
270{
271 vector<string> v;
272 string pushstring = "";
273 int index;
274 float writestring;
275 for ( int i = 0; i < h_numSFperElem( attype ); ++i )
276 {
277 index = SF( attype, i, 13 );
278 // TODO: improve function
279 for ( int j = 1; j < 12; ++j )
280 {
281 writestring = SF( attype, index, j );
282 pushstring += to_string( writestring ) + " ";
283 }
284 pushstring += "\n";
285 }
286 v.push_back( pushstring );
287
288 return v;
289}
290
291template <class t_SF, class t_SFscaling, class h_t_int>
292vector<string>
293ElementCabana::infoSymmetryFunctionScaling( ScalingType scalingType, t_SF SF,
294 t_SFscaling SFscaling, int attype,
295 h_t_int h_numSFperElem ) const
296{
297 vector<string> v;
298 int index;
299 for ( int k = 0; k < h_numSFperElem( attype ); ++k )
300 {
301 index = SF( attype, k, 13 );
302 v.push_back( scalingLine( scalingType, SFscaling, attype, index ) );
303 }
304 return v;
305}
306
307template <class t_SF, class t_SFGmemberlist, class h_t_int>
308void ElementCabana::setupSymmetryFunctionGroups( t_SF SF,
309 t_SFGmemberlist SFGmemberlist,
310 int attype, h_t_int h_numSFperElem,
311 h_t_int h_numSFGperElem,
312 int maxSFperElem )
313{
314 int num_group = h_numSFperElem.extent( 0 );
315 h_t_int h_numGR( "RadialCounter", num_group );
316 h_t_int h_numGA( "AngularCounter", num_group );
317 int SFindex;
318 for ( int k = 0; k < h_numSFperElem( attype ); ++k )
319 {
320 bool createNewGroup = true;
321 SFindex = SF( attype, k, 13 );
322 for ( int l = 0; l < h_numSFGperElem( attype ); ++l )
323 {
324 if ( ( SF( attype, SFindex, 0 ) ==
325 SF( attype, SFGmemberlist( attype, l, 0 ),
326 0 ) ) && // same ec
327 ( SF( attype, SFindex, 2 ) ==
328 SF( attype, SFGmemberlist( attype, l, 0 ),
329 2 ) ) && // same e1
330 ( SF( attype, SFindex, 3 ) ==
331 SF( attype, SFGmemberlist( attype, l, 0 ),
332 3 ) ) && // same e2
333 ( SF( attype, SFindex, 7 ) ==
334 SF( attype, SFGmemberlist( attype, l, 0 ),
335 7 ) ) && // same rc
336 ( SF( attype, SFindex, 11 ) ==
337 SF( attype, SFGmemberlist( attype, l, 0 ),
338 11 ) ) && // same cutoffType
339 ( SF( attype, SFindex, 12 ) ==
340 SF( attype, SFGmemberlist( attype, l, 0 ),
341 12 ) ) ) // same cutoffAlpha
342 {
343 createNewGroup = false;
344 if ( SF( attype, SFindex, 1 ) == 2 )
345 {
346 SFGmemberlist( attype, l, h_numGR( l ) ) = SFindex;
347 h_numGR( l )++;
348 SFGmemberlist( attype, l, maxSFperElem )++;
349 break;
350 }
351
352 else if ( SF( attype, SFindex, 1 ) == 3 )
353 {
354 SFGmemberlist( attype, l, h_numGA( l ) ) = SFindex;
355 h_numGA( l )++;
356 SFGmemberlist( attype, l, maxSFperElem )++;
357 break;
358 }
359 }
360 }
361
362 if ( createNewGroup )
363 {
364 int l = h_numSFGperElem( attype );
365 h_numSFGperElem( attype )++;
366 if ( SF( attype, SFindex, 1 ) == 2 )
367 {
368 SFGmemberlist( attype, l, 0 ) = SFindex;
369 if ( l >= (int)h_numGR.extent( 0 ) )
370 Kokkos::resize( h_numGR, l+1 );
371 h_numGR( l ) = 1;
372 SFGmemberlist( attype, l, maxSFperElem )++;
373 }
374 else if ( SF( attype, SFindex, 1 ) == 3 )
375 {
376 SFGmemberlist( attype, l, 0 ) = SFindex;
377 if ( l >= (int)h_numGA.extent( 0 ) )
378 Kokkos::resize( h_numGA, l+1 );
379 h_numGA( l ) = 1;
380 SFGmemberlist( attype, l, maxSFperElem )++;
381 }
382 }
383 }
384
385 return;
386}
387
388template <class t_SF, class t_SFGmemberlist, class h_t_int>
389vector<string>
390ElementCabana::infoSymmetryFunctionGroups( t_SF SF, t_SFGmemberlist SFGmemberlist,
391 int attype, h_t_int h_numSFGperElem ) const
392{
393 vector<string> v;
394 string pushstring = "";
395 for ( int groupIndex = 0; groupIndex < h_numSFGperElem( attype );
396 ++groupIndex )
397 {
398 // TODO: improve function
399 for ( int j = 0; j < 8; ++j )
400 pushstring +=
401 to_string(
402 SF( attype, SFGmemberlist( attype, groupIndex, 0 ), j ) ) +
403 " ";
404 pushstring += "\n";
405 }
406 v.push_back( pushstring );
407
408 return v;
409}
410
411template <class t_SF, class h_t_int>
412void ElementCabana::setCutoffFunction( CutoffFunction::CutoffType const cutoffType,
413 double const cutoffAlpha, t_SF SF, int attype,
414 h_t_int h_numSFperElem )
415{
416 for ( int k = 0; k < h_numSFperElem( attype ); ++k )
417 {
418 SF( attype, k, 10 ) = cutoffType;
419 SF( attype, k, 11 ) = cutoffAlpha;
420 }
421 return;
422}
423
424template <class t_SF, class t_SFscaling, class h_t_int>
425void ElementCabana::setScaling( ScalingType scalingType,
426 vector<string> const &statisticsLine, double Smin,
427 double Smax, t_SF SF, t_SFscaling SFscaling,
428 int attype, h_t_int h_numSFperElem ) const
429{
430 int index;
431 for ( int k = 0; k < h_numSFperElem( attype ); ++k )
432 {
433 index = SF( attype, k, 13 );
434 setScalingType( scalingType, statisticsLine.at( k ), Smin, Smax,
435 SFscaling, attype, index );
436 }
437 // TODO: groups
438 // for (int k = 0; k < h_numSFperElem(attype); ++k)
439 // setScalingFactors(SF,attype,k);
440
441 return;
442}
443
444template <class t_SF>
445size_t ElementCabana::getMinNeighbors( int attype, t_SF SF, int nSF ) const
446{
447 // get max number of minNeighbors
448 size_t global_minNeighbors = 0;
449 size_t minNeighbors = 0;
450 int SFtype;
451 for ( int k = 0; k < nSF; ++k )
452 {
453 SFtype = SF( attype, k, 1 );
454 if ( SFtype == 2 )
455 minNeighbors = 1;
456 else if ( SFtype == 3 )
457 minNeighbors = 2;
458 global_minNeighbors = max( minNeighbors, global_minNeighbors );
459 }
460
461 return global_minNeighbors;
462}
463
464template <class t_SF, class h_t_int>
465double ElementCabana::getMinCutoffRadius( t_SF SF, int attype,
466 h_t_int h_numSFperElem ) const
467{
468 double minCutoffRadius = numeric_limits<double>::max();
469
470 for ( int k = 0; k < h_numSFperElem( attype ); ++k )
471 minCutoffRadius = min( SF( attype, k, 7 ), minCutoffRadius );
472
473 return minCutoffRadius;
474}
475
476template <class t_SF, class h_t_int>
477double ElementCabana::getMaxCutoffRadius( t_SF SF, int attype,
478 h_t_int h_numSFperElem ) const
479{
480 double maxCutoffRadius = 0.0;
481
482 for ( int k = 0; k < h_numSFperElem( attype ); ++k )
483 maxCutoffRadius = max( SF( attype, k, 7 ), maxCutoffRadius );
484
485 return maxCutoffRadius;
486}
487
488// TODO:add functionality
489/*
490void ElementCabana::updateSymmetryFunctionStatistics
491*/
492
493}
nnp::CutoffFunction::CutoffType
CutoffType
List of available cutoff function types.
Definition: CutoffFunction.h:43
nnp::ElementCabana::setScaling
void setScaling(ScalingType scalingType, std::vector< std::string > const &statisticsLine, double minS, double maxS, t_SF SF, t_SFscaling SFscaling, int attype, h_t_int h_numSFperElem) const
Set scaling of all symmetry functions.
Definition: ElementCabana_impl.h:425
nnp::ElementCabana::~ElementCabana
~ElementCabana()
Destructor.
Definition: ElementCabana_impl.h:39
nnp::ElementCabana::setScalingType
void setScalingType(ScalingType scalingType, std::string statisticsLine, double Smin, double Smax, t_SFscaling SFscaling, int attype, int k) const
Symmetry function statistics.
Definition: ElementCabana.h:286
nnp::ElementCabana::index
std::size_t index
Global index of this element.
Definition: Element.h:255
nnp::ElementCabana::setCutoffFunction
void setCutoffFunction(CutoffFunction::CutoffType const cutoffType, double const cutoffAlpha, t_SF SF, int attype, h_t_int h_numSFperElem)
Set cutoff function for all symmetry functions.
Definition: ElementCabana_impl.h:412
nnp::ElementCabana::compareSF
bool compareSF(t_SF SF, int attype, int index1, int index2)
Print symmetry function parameter value information.
Definition: ElementCabana_impl.h:205
nnp::ElementCabana::ElementCabana
ElementCabana(std::size_t const index)
Constructor using index.
Definition: ElementCabana_impl.h:31
nnp::ElementCabana::atomicEnergyOffset
double atomicEnergyOffset
Offset energy for every atom of this element.
Definition: Element.h:259
nnp::ElementCabana::scalingLine
std::string scalingLine(ScalingType scalingType, t_SFscaling SFscaling, int attype, int k) const
Print scaling for one symmetry function.
Definition: ElementCabana.h:327
nnp::ElementCabana::addSymmetryFunction
void addSymmetryFunction(std::string const &parameters, std::vector< std::string > elementStrings, int attype, t_SF SF, double convLength, h_t_int h_numSFperElem)
Add one symmetry function.
Definition: ElementCabana_impl.h:42
nnp::Element
Contains element-specific data.
Definition: Element.h:39
nnp::Element::infoSymmetryFunctionParameters
std::vector< std::string > infoSymmetryFunctionParameters() const
Print symmetry function parameter value information.
Definition: Element.cpp:168
nnp::Element::sortSymmetryFunctions
void sortSymmetryFunctions()
Sort all symmetry function.
Definition: Element.cpp:154
nnp::Element::infoSymmetryFunctionGroups
std::vector< std::string > infoSymmetryFunctionGroups() const
Print symmetry function group info.
Definition: Element.cpp:313
nnp::Element::infoSymmetryFunctionScaling
std::vector< std::string > infoSymmetryFunctionScaling() const
Print symmetry function scaling information.
Definition: Element.cpp:181
nnp::Element::getMinCutoffRadius
double getMinCutoffRadius() const
Get minimum cutoff radius of all symmetry functions.
Definition: Element.cpp:397
nnp::Element::setupSymmetryFunctionGroups
void setupSymmetryFunctionGroups()
Set up symmetry function groups.
Definition: Element.cpp:194
nnp::Element::getMinNeighbors
std::size_t getMinNeighbors() const
Get maximum of required minimum number of neighbors for all symmetry functions for this element.
Definition: Element.cpp:384
nnp::Element::getMaxCutoffRadius
double getMaxCutoffRadius() const
Get maximum cutoff radius of all symmetry functions.
Definition: Element.cpp:414
nnp
Definition: Atom.h:29
nnp::split
vector< string > split(string const &input, char delimiter)
Split string at each delimiter.
Definition: utility.cpp:33
nnp::reduce
string reduce(string const &line, string const &whitespace, string const &fill)
Replace multiple whitespaces with fill.
Definition: utility.cpp:60
ScalingType
ScalingType
Definition: typesCabana.h:26