1 |
gezelter |
1862 |
/* |
2 |
|
|
* Copyright (c) 2013 The University of Notre Dame. All Rights Reserved. |
3 |
|
|
* |
4 |
|
|
* The University of Notre Dame grants you ("Licensee") a |
5 |
|
|
* non-exclusive, royalty free, license to use, modify and |
6 |
|
|
* redistribute this software in source and binary code form, provided |
7 |
|
|
* that the following conditions are met: |
8 |
|
|
* |
9 |
|
|
* 1. Redistributions of source code must retain the above copyright |
10 |
|
|
* notice, this list of conditions and the following disclaimer. |
11 |
|
|
* |
12 |
|
|
* 2. Redistributions in binary form must reproduce the above copyright |
13 |
|
|
* notice, this list of conditions and the following disclaimer in the |
14 |
|
|
* documentation and/or other materials provided with the |
15 |
|
|
* distribution. |
16 |
|
|
* |
17 |
|
|
* This software is provided "AS IS," without a warranty of any |
18 |
|
|
* kind. All express or implied conditions, representations and |
19 |
|
|
* warranties, including any implied warranty of merchantability, |
20 |
|
|
* fitness for a particular purpose or non-infringement, are hereby |
21 |
|
|
* excluded. The University of Notre Dame and its licensors shall not |
22 |
|
|
* be liable for any damages suffered by licensee as a result of |
23 |
|
|
* using, modifying or distributing the software or its |
24 |
|
|
* derivatives. In no event will the University of Notre Dame or its |
25 |
|
|
* licensors be liable for any lost revenue, profit or data, or for |
26 |
|
|
* direct, indirect, special, consequential, incidental or punitive |
27 |
|
|
* damages, however caused and regardless of the theory of liability, |
28 |
|
|
* arising out of the use of or inability to use software, even if the |
29 |
|
|
* University of Notre Dame has been advised of the possibility of |
30 |
|
|
* such damages. |
31 |
|
|
* |
32 |
|
|
* SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your |
33 |
|
|
* research, please cite the appropriate papers when you publish your |
34 |
|
|
* work. Good starting points are: |
35 |
|
|
* |
36 |
|
|
* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). |
37 |
|
|
* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). |
38 |
|
|
* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). |
39 |
|
|
* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). |
40 |
|
|
* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). |
41 |
|
|
*/ |
42 |
|
|
|
43 |
|
|
#include "clusters/Icosahedron.hpp" |
44 |
|
|
|
45 |
|
|
using namespace std; |
46 |
|
|
|
47 |
|
|
namespace OpenMD { |
48 |
|
|
|
49 |
|
|
Icosahedron::Icosahedron() { |
50 |
|
|
Basis.clear(); |
51 |
|
|
Edges.clear(); |
52 |
|
|
Facets.clear(); |
53 |
|
|
Points.clear(); |
54 |
|
|
|
55 |
|
|
// |
56 |
|
|
// Initialize Basis vectors. |
57 |
|
|
// |
58 |
|
|
const RealType HGR = ( sqrt(5.0) + 1.0 ) / 4.0; // half of the golden ratio |
59 |
|
|
|
60 |
|
|
Basis.push_back( Vector3d( HGR, 0.0, 0.5 )); |
61 |
|
|
Basis.push_back( Vector3d( HGR, 0.0, -0.5 )); |
62 |
|
|
Basis.push_back( Vector3d( 0.5, HGR, 0.0 )); |
63 |
|
|
Basis.push_back( Vector3d( -0.5, HGR, 0.0 )); |
64 |
|
|
Basis.push_back( Vector3d( 0.0, 0.5, HGR )); |
65 |
|
|
Basis.push_back( Vector3d( 0.0, -0.5, HGR )); |
66 |
|
|
Basis.push_back( Vector3d( 0.5, -HGR, 0.0 )); |
67 |
|
|
Basis.push_back( Vector3d( 0.0, 0.5, -HGR )); |
68 |
|
|
Basis.push_back( Vector3d( -HGR, 0.0, 0.5 )); |
69 |
|
|
Basis.push_back( Vector3d( 0.0, -0.5, -HGR )); |
70 |
|
|
Basis.push_back( Vector3d( -HGR, 0.0, -0.5 )); |
71 |
|
|
Basis.push_back( Vector3d( -0.5, -HGR, 0.0 )); |
72 |
|
|
|
73 |
|
|
// |
74 |
|
|
// Initialize 30 edges |
75 |
|
|
// |
76 |
|
|
|
77 |
|
|
Edges.push_back(std::make_pair(0, 1)); |
78 |
|
|
Edges.push_back(std::make_pair(0, 2)); |
79 |
|
|
Edges.push_back(std::make_pair(0, 4)); |
80 |
|
|
Edges.push_back(std::make_pair(0, 5)); |
81 |
|
|
Edges.push_back(std::make_pair(0, 6)); |
82 |
|
|
|
83 |
|
|
Edges.push_back(std::make_pair(10, 3)); |
84 |
|
|
Edges.push_back(std::make_pair(10, 7)); |
85 |
|
|
Edges.push_back(std::make_pair(10, 8)); |
86 |
|
|
Edges.push_back(std::make_pair(10, 9)); |
87 |
|
|
Edges.push_back(std::make_pair(10, 11)); |
88 |
|
|
|
89 |
|
|
Edges.push_back(std::make_pair(1, 2)); |
90 |
|
|
Edges.push_back(std::make_pair(1, 6)); |
91 |
|
|
Edges.push_back(std::make_pair(1, 7)); |
92 |
|
|
Edges.push_back(std::make_pair(1, 9)); |
93 |
|
|
|
94 |
|
|
Edges.push_back(std::make_pair(8, 3)); |
95 |
|
|
Edges.push_back(std::make_pair(8, 4)); |
96 |
|
|
Edges.push_back(std::make_pair(8, 5)); |
97 |
|
|
Edges.push_back(std::make_pair(8, 11)); |
98 |
|
|
|
99 |
|
|
Edges.push_back(std::make_pair(2, 3)); |
100 |
|
|
Edges.push_back(std::make_pair(2, 4)); |
101 |
|
|
Edges.push_back(std::make_pair(2, 7)); |
102 |
|
|
|
103 |
|
|
Edges.push_back(std::make_pair(11, 5)); |
104 |
|
|
Edges.push_back(std::make_pair(11, 6)); |
105 |
|
|
Edges.push_back(std::make_pair(11, 9)); |
106 |
|
|
|
107 |
|
|
Edges.push_back(std::make_pair(6, 5)); |
108 |
|
|
Edges.push_back(std::make_pair(6, 9)); |
109 |
|
|
|
110 |
|
|
Edges.push_back(std::make_pair(3, 4)); |
111 |
|
|
Edges.push_back(std::make_pair(3, 7)); |
112 |
|
|
|
113 |
|
|
Edges.push_back(std::make_pair(7, 9)); |
114 |
|
|
|
115 |
|
|
Edges.push_back(std::make_pair(5, 4)); |
116 |
|
|
|
117 |
|
|
// |
118 |
|
|
// Initialize 20 facets |
119 |
|
|
// |
120 |
|
|
|
121 |
|
|
Facets.push_back(make_tuple3(0, 1, 2)); |
122 |
|
|
Facets.push_back(make_tuple3(0, 2, 4)); |
123 |
|
|
Facets.push_back(make_tuple3(0, 4, 5)); |
124 |
|
|
Facets.push_back(make_tuple3(0, 5, 6)); |
125 |
|
|
Facets.push_back(make_tuple3(0, 1, 6)); |
126 |
|
|
|
127 |
|
|
Facets.push_back(make_tuple3(10, 3, 7)); |
128 |
|
|
Facets.push_back(make_tuple3(10, 3, 8)); |
129 |
|
|
Facets.push_back(make_tuple3(10, 8, 11)); |
130 |
|
|
Facets.push_back(make_tuple3(10, 9, 11)); |
131 |
|
|
Facets.push_back(make_tuple3(10, 7, 9)); |
132 |
|
|
|
133 |
|
|
Facets.push_back(make_tuple3(1, 2, 7)); |
134 |
|
|
Facets.push_back(make_tuple3(1, 7, 9)); |
135 |
|
|
Facets.push_back(make_tuple3(1, 6, 9)); |
136 |
|
|
|
137 |
|
|
Facets.push_back(make_tuple3(8, 5, 11)); |
138 |
|
|
Facets.push_back(make_tuple3(8, 4, 5)); |
139 |
|
|
Facets.push_back(make_tuple3(8, 3, 4)); |
140 |
|
|
|
141 |
|
|
Facets.push_back(make_tuple3(2, 3, 7)); |
142 |
|
|
Facets.push_back(make_tuple3(2, 3, 4)); |
143 |
|
|
|
144 |
|
|
Facets.push_back(make_tuple3(11, 5, 6)); |
145 |
|
|
Facets.push_back(make_tuple3(11, 6, 9)); |
146 |
|
|
} |
147 |
|
|
|
148 |
|
|
Icosahedron::~Icosahedron() { |
149 |
|
|
Facets.clear(); |
150 |
|
|
Edges.clear(); |
151 |
|
|
Basis.clear(); |
152 |
|
|
Points.clear(); |
153 |
|
|
} |
154 |
|
|
|
155 |
|
|
int Icosahedron::getNpoints( int n ) { |
156 |
|
|
int count=0; |
157 |
|
|
for( int i = 0; i <= n; i++ ) count += np( i ); |
158 |
|
|
return count; |
159 |
|
|
} |
160 |
|
|
|
161 |
|
|
int Icosahedron::np( int n ) { |
162 |
|
|
if( n<0 ) return -1; |
163 |
|
|
else if( n==0 ) return 1; |
164 |
|
|
else if( n==1 ) return 12; |
165 |
|
|
else if( n==2 ) return 42; |
166 |
|
|
else { |
167 |
|
|
int count = 0; |
168 |
|
|
count += 12; // edge particles |
169 |
|
|
count += (n-1)*30; // side particles |
170 |
|
|
for( int i = 1; i <= n-2; i++ ) count += i*20; // body particles |
171 |
|
|
return count; |
172 |
|
|
} |
173 |
|
|
} |
174 |
|
|
|
175 |
|
|
vector<Vector3d> Icosahedron::ih( int n ) { |
176 |
|
|
|
177 |
|
|
if( n < 0 ) return Points; |
178 |
|
|
|
179 |
|
|
if( n==0 ) { |
180 |
|
|
|
181 |
|
|
// center particle only |
182 |
|
|
|
183 |
|
|
Points.push_back(Vector3d( 0.0, 0.0, 0.0 )); |
184 |
|
|
return Points; |
185 |
|
|
} |
186 |
|
|
|
187 |
|
|
// |
188 |
|
|
// Generate edge particles |
189 |
|
|
// |
190 |
|
|
|
191 |
|
|
for( vector<Vector3d>::iterator i = Basis.begin(); i != Basis.end(); ++i ) { |
192 |
|
|
|
193 |
|
|
Points.push_back( (*i) * RealType(n) ); |
194 |
|
|
} |
195 |
|
|
|
196 |
|
|
// |
197 |
|
|
// Generate side particles |
198 |
|
|
// |
199 |
|
|
|
200 |
|
|
if( n<2 ) return Points; |
201 |
|
|
|
202 |
|
|
for( vector<pair<int,int> >::iterator i=Edges.begin(); |
203 |
|
|
i != Edges.end(); ++i ) { |
204 |
|
|
|
205 |
|
|
Vector3d e1 = Basis[ (*i).first ] * RealType(n); |
206 |
|
|
Vector3d e2 = Basis[ (*i).second ] * RealType(n); |
207 |
|
|
|
208 |
|
|
for( int j = 1; j <= n-1; j++ ) { |
209 |
|
|
Points.push_back( e1 + (e2-e1) * RealType(j) / RealType(n)); |
210 |
|
|
} |
211 |
|
|
} |
212 |
|
|
|
213 |
|
|
// |
214 |
|
|
// Generate body particles |
215 |
|
|
// |
216 |
|
|
|
217 |
|
|
if( n<3 ) return Points; |
218 |
|
|
|
219 |
|
|
for( vector<tuple3<int,int,int> >::iterator i = Facets.begin(); |
220 |
|
|
i != Facets.end(); ++i) { |
221 |
|
|
|
222 |
|
|
Vector3d e1 = Basis[ (*i).first ] * RealType(n); |
223 |
|
|
Vector3d e2 = Basis[ (*i).second ] * RealType(n); |
224 |
|
|
Vector3d e3 = Basis[ (*i).third ] * RealType(n); |
225 |
|
|
|
226 |
|
|
for( int j=1; j<=n-2; j++ ) { |
227 |
|
|
|
228 |
|
|
Vector3d v1 = e1 + (e2-e1) * RealType(j+1) / RealType(n); |
229 |
|
|
Vector3d v2 = e1 + (e3-e1) * RealType(j+1) / RealType(n); |
230 |
|
|
|
231 |
|
|
for( int k=1; k<=j; k++ ) { |
232 |
|
|
Points.push_back(v1 + (v2-v1) * RealType(k) / RealType(j+1)); |
233 |
|
|
} |
234 |
|
|
} |
235 |
|
|
} |
236 |
|
|
return Points; |
237 |
|
|
} |
238 |
|
|
|
239 |
|
|
vector<Vector3d> Icosahedron::getPoints(int nshells) { |
240 |
|
|
//generate the coordinates |
241 |
|
|
for( int i = 0; i <= nshells; i++ ) ih( i ); |
242 |
|
|
return Points; |
243 |
|
|
} |
244 |
|
|
} |