1 // -*- C++ -*-
2 // KimAtoms.cpp: Interface to KIM pretending to be the atoms.
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 // This is a reimplementation of the ASAP Atoms object. It does not
31 // inherit from the usual Atoms object, as Atoms.h include Python.h
32 // which is not available in OpenKIM.
33
34 #include "KimAtoms.h"
35 #include "Exception.h"
36 #include "Debug.h"
37 #include <math.h>
38
KimAtoms()39 KimAtoms::KimAtoms()
40 {
41 CONSTRUCTOR;
42 refcount = 1;
43 counter = 2;
44 count_inverse_cell = 1;
45 }
46
~KimAtoms()47 KimAtoms::~KimAtoms()
48 {
49 DESTRUCTOR;
50 }
51
ReInit(KIM::ModelComputeArguments const * const modelComputeArguments,int nAtoms,double const * pos,int const * z,int const * contributes)52 void KimAtoms::ReInit(KIM::ModelComputeArguments const * const modelComputeArguments,
53 int nAtoms,
54 double const *pos,
55 int const *z,
56 int const *contributes)
57 {
58 this->modelComputeArguments = modelComputeArguments;
59 this->nAtoms = nAtoms;
60 positions.resize(nAtoms);
61 numbers.resize(nAtoms);
62 this->contributes = contributes;
63 for (int i = 0; i < nAtoms; i++)
64 {
65 positions[i] = ((Vec *)pos)[i];
66 numbers[i] = (asap_z_int) z[i];
67 }
68 counter++;
69 for (int i = 0; i < 3; i++)
70 {
71 for (int j = 0; j < 3; j++)
72 cell[i][j] = (i == j) * 50.0; // Fake number as it is not really used.
73 }
74 }
75
GetListOfElements(set<int> & elements) const76 void KimAtoms::GetListOfElements(set<int> &elements) const
77 {
78 DEBUGPRINT;
79 const asap_z_int *atomicnumbers = GetAtomicNumbers();
80
81 elements.clear();
82 for (int i = 0; i < nAtoms; i++)
83 {
84 int z = atomicnumbers[i];
85 if (elements.find(z) == elements.end())
86 elements.insert(z);
87 }
88 DEBUGPRINT;
89 }
90
91 /// Get the cartesian momenta
GetMomenta()92 const Vec *KimAtoms::GetMomenta()
93 {
94 throw AsapError("GetMomenta called in KIM mode.");
95 return NULL;
96 }
97
98
99 /// Get the cartesian momenta
GetMasses()100 const double *KimAtoms::GetMasses()
101 {
102 throw AsapError("GetMasses called in KIM mode.");
103 return NULL;
104 }
105
GetVolume() const106 double KimAtoms::GetVolume() const
107 {
108 DEBUGPRINT;
109 double det;
110 det = -cell[0][2]*cell[1][1]*cell[2][0] +
111 cell[0][1]*cell[1][2]*cell[2][0] +
112 cell[0][2]*cell[1][0]*cell[2][1] -
113 cell[0][0]*cell[1][2]*cell[2][1] -
114 cell[0][1]*cell[1][0]*cell[2][2] +
115 cell[0][0]*cell[1][1]*cell[2][2];
116 DEBUGPRINT;
117 return fabs(det);
118 }
119
GetPositions(vector<Vec> & pos,bool ghosts) const120 void KimAtoms::GetPositions(vector<Vec> &pos, bool ghosts /* = false */) const
121 {
122 DEBUGPRINT;
123 pos.clear();
124 int n = nAtoms;
125 pos.reserve(n + n/25); // 4% extra
126 for (int i = 0; i < n; i++)
127 pos[i] = positions[i];
128
129 DEBUGPRINT;
130 }
131
GetScaledPositions(vector<Vec> & pos,bool ghosts)132 void KimAtoms::GetScaledPositions(vector<Vec> &pos, bool ghosts /* = false */)
133 {
134 DEBUGPRINT;
135 int n = nAtoms;
136 const Vec *inv = GetInverseCell();
137 if (pos.capacity() < n)
138 pos.reserve(n + n/25); // Reserve 4% extra.
139 pos.resize(n);
140 for (int i = 0; i < n; i++)
141 for (int j = 0; j < 3; j++)
142 pos[i][j] = positions[i][0] * inv[0][j]
143 + positions[i][1] * inv[1][j]
144 + positions[i][2] * inv[2][j];
145 DEBUGPRINT;
146 }
147
GetScaledPositions(vector<Vec> & scaledpos,const set<int> & which)148 void KimAtoms::GetScaledPositions(vector<Vec> &scaledpos, const set<int> &which)
149 {
150 DEBUGPRINT;
151 assert(scaledpos.size() == which.size());
152 const Vec *inv = GetInverseCell();
153 vector<Vec>::iterator spi = scaledpos.begin();
154 for (set<int>::const_iterator i = which.begin(); i != which.end(); ++i,++spi)
155 for (int j = 0; j < 3; j++)
156 (*spi)[j] = positions[*i][0] * inv[0][j]
157 + positions[*i][1] * inv[1][j]
158 + positions[*i][2] * inv[2][j];
159 DEBUGPRINT;
160 }
161
GetCellHeights()162 const double *KimAtoms::GetCellHeights()
163 {
164 DEBUGPRINT;
165 if (count_inverse_cell < counter)
166 invert_cell();
167 return heights;
168 }
169
GetInverseCell()170 const Vec *KimAtoms::GetInverseCell()
171 {
172 DEBUGPRINT;
173 if (count_inverse_cell < counter)
174 invert_cell();
175 return inverse;
176 }
177
invert_cell()178 void KimAtoms::invert_cell()
179 {
180 DEBUGPRINT;
181 count_inverse_cell = counter;
182 double determinant = Cross(cell[0], cell[1]) * cell[2];
183 // Find heights
184 for (int i = 0; i < 3; i++)
185 {
186 Vec inv = Cross(cell[(i + 1) % 3], cell[(i + 2) % 3]);
187 heights[i] = fabs(determinant) / sqrt(Length2(inv));
188 }
189 // Invert matrix. I_i,j = { C_j-1,i-1 C_j+1,i+1 - C_j+1,i-1 C_j-1,i+1 } / det
190 for (int i = 0; i < 3; i++)
191 {
192 int ip = (i + 1) % 3;
193 int im = (i + 2) % 3;
194 for (int j = 0; j < 3; j++)
195 {
196 int jp = (j + 1) % 3;
197 int jm = (j + 2) % 3;
198 inverse[i][j] = (cell[jm][im]*cell[jp][ip] -
199 cell[jp][im]*cell[jm][ip]) / determinant;
200 }
201 }
202 DEBUGPRINT;
203 }
204
SetDiagonalCell(double d[3])205 void KimAtoms::SetDiagonalCell(double d[3])
206 {
207 for (int i = 0; i < 3; i++)
208 for (int j = 0; j < 3; j++)
209 if (i == j)
210 cell[i][j] = d[i];
211 else
212 cell[i][j] = 0.0;
213 count_inverse_cell = 0; // Must be recalculated
214 }
215
216