src/clusters/Icosahedron.cpp

/*
 * Copyright (c) 2013 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 *
 * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
 * research, please cite the appropriate papers when you publish your
 * work.  Good starting points are:
 *                                                                      
 * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
 * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
 * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
 * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
 * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
 */

#include "clusters/Icosahedron.hpp"

using namespace std;

namespace OpenMD {

  Icosahedron::Icosahedron() {
    Basis.clear();
    Edges.clear();
    Facets.clear();
    Points.clear();
    
    //
    // Initialize Basis vectors.
    //
    const RealType HGR = ( sqrt(5.0) + 1.0 ) / 4.0; // half of the golden ratio
    
    Basis.push_back( Vector3d(  HGR,  0.0,  0.5 ));
    Basis.push_back( Vector3d(  HGR,  0.0, -0.5 ));
    Basis.push_back( Vector3d(  0.5,  HGR,  0.0 ));
    Basis.push_back( Vector3d( -0.5,  HGR,  0.0 ));
    Basis.push_back( Vector3d(  0.0,  0.5,  HGR ));
    Basis.push_back( Vector3d(  0.0, -0.5,  HGR ));
    Basis.push_back( Vector3d(  0.5, -HGR,  0.0 ));
    Basis.push_back( Vector3d(  0.0,  0.5, -HGR ));
    Basis.push_back( Vector3d( -HGR,  0.0,  0.5 ));
    Basis.push_back( Vector3d(  0.0, -0.5, -HGR ));
    Basis.push_back( Vector3d( -HGR,  0.0, -0.5 ));
    Basis.push_back( Vector3d( -0.5, -HGR,  0.0 ));
    
    //
    // Initialize 30 edges
    //
    
    Edges.push_back(std::make_pair(0, 1));
    Edges.push_back(std::make_pair(0, 2));
    Edges.push_back(std::make_pair(0, 4));
    Edges.push_back(std::make_pair(0, 5));
    Edges.push_back(std::make_pair(0, 6));
    
    Edges.push_back(std::make_pair(10, 3));
    Edges.push_back(std::make_pair(10, 7));
    Edges.push_back(std::make_pair(10, 8));
    Edges.push_back(std::make_pair(10, 9));
    Edges.push_back(std::make_pair(10, 11));
  
    Edges.push_back(std::make_pair(1, 2));
    Edges.push_back(std::make_pair(1, 6));
    Edges.push_back(std::make_pair(1, 7));
    Edges.push_back(std::make_pair(1, 9));
    
    Edges.push_back(std::make_pair(8, 3));
    Edges.push_back(std::make_pair(8, 4));
    Edges.push_back(std::make_pair(8, 5));
    Edges.push_back(std::make_pair(8, 11));
    
    Edges.push_back(std::make_pair(2, 3));
    Edges.push_back(std::make_pair(2, 4));
    Edges.push_back(std::make_pair(2, 7));
    
    Edges.push_back(std::make_pair(11, 5));
    Edges.push_back(std::make_pair(11, 6));
    Edges.push_back(std::make_pair(11, 9));
    
    Edges.push_back(std::make_pair(6, 5));
    Edges.push_back(std::make_pair(6, 9));
    
    Edges.push_back(std::make_pair(3, 4));
    Edges.push_back(std::make_pair(3, 7));
    
    Edges.push_back(std::make_pair(7, 9));
    
    Edges.push_back(std::make_pair(5, 4));
    
    //
    // Initialize 20 facets
    //
    
    Facets.push_back(make_tuple3(0, 1, 2));
    Facets.push_back(make_tuple3(0, 2, 4));
    Facets.push_back(make_tuple3(0, 4, 5));
    Facets.push_back(make_tuple3(0, 5, 6));
    Facets.push_back(make_tuple3(0, 1, 6));
    
    Facets.push_back(make_tuple3(10, 3, 7));
    Facets.push_back(make_tuple3(10, 3, 8));
    Facets.push_back(make_tuple3(10, 8, 11));
    Facets.push_back(make_tuple3(10, 9, 11));
    Facets.push_back(make_tuple3(10, 7, 9));
    
    Facets.push_back(make_tuple3(1, 2, 7));
    Facets.push_back(make_tuple3(1, 7, 9));
    Facets.push_back(make_tuple3(1, 6, 9));
    
    Facets.push_back(make_tuple3(8, 5, 11));
    Facets.push_back(make_tuple3(8, 4, 5));
    Facets.push_back(make_tuple3(8, 3, 4));
    
    Facets.push_back(make_tuple3(2, 3, 7));
    Facets.push_back(make_tuple3(2, 3, 4));
    
    Facets.push_back(make_tuple3(11, 5, 6));
    Facets.push_back(make_tuple3(11, 6, 9));
  }
  
  Icosahedron::~Icosahedron() {
    Facets.clear();
    Edges.clear();
    Basis.clear();
    Points.clear();
  }
  
  int Icosahedron::getNpoints( int n ) {
    int count=0;
    for( int i = 0; i <= n; i++ ) count += np( i );
    return count;
  }

  int Icosahedron::np( int n ) {
    if( n<0 ) return -1;
    else if( n==0 ) return 1;
    else if( n==1 ) return 12;
    else if( n==2 ) return 42;
    else {
      int count = 0;    
      count += 12;   // edge particles
      count += (n-1)*30;   // side particles
      for( int i = 1; i <= n-2; i++ ) count += i*20;   // body particles
      return count;
    }
  }

  vector<Vector3d> Icosahedron::ih( int n ) {
    
    if( n < 0 ) return Points; 
    
    if( n==0 ) {

      // center particle only

      Points.push_back(Vector3d( 0.0, 0.0, 0.0 ));
      return Points;
    } 
    
    //
    // Generate edge particles
    //

    for( vector<Vector3d>::iterator i = Basis.begin(); i != Basis.end(); ++i ) {
      
      Points.push_back( (*i) * RealType(n) );
    }
    
    //
    // Generate side particles
    //
    
    if( n<2 ) return Points;
    
    for( vector<pair<int,int> >::iterator i=Edges.begin(); 
         i != Edges.end(); ++i ) {
      
      Vector3d e1 = Basis[ (*i).first  ] * RealType(n);
      Vector3d e2 = Basis[ (*i).second ] * RealType(n);
      
      for( int j = 1; j <= n-1; j++ ) {
        Points.push_back( e1 + (e2-e1) * RealType(j) / RealType(n));
      }      
    }
    
    //
    // Generate body particles
    //
    
    if( n<3 ) return Points;
    
    for( vector<tuple3<int,int,int> >::iterator i = Facets.begin();
         i != Facets.end(); ++i) {
      
      Vector3d e1 = Basis[ (*i).first  ] * RealType(n);
      Vector3d e2 = Basis[ (*i).second ] * RealType(n);
      Vector3d e3 = Basis[ (*i).third  ] * RealType(n);
      
      for( int j=1; j<=n-2; j++ ) {
        
        Vector3d v1 = e1 + (e2-e1) * RealType(j+1) / RealType(n);
        Vector3d v2 = e1 + (e3-e1) * RealType(j+1) / RealType(n);
        
        for( int k=1; k<=j; k++ ) {
          Points.push_back(v1 + (v2-v1) * RealType(k) / RealType(j+1));
        }
      }
    }
    return Points;
  }
  
  vector<Vector3d> Icosahedron::getPoints(int nshells) {
    //generate the coordinates
    for( int i = 0; i <= nshells; i++ ) ih( i );
    return Points;
  }
}
Revision:	1879
Committed:	Sun Jun 16 15:15:42 2013 UTC (12 years, 4 months ago) by gezelter
File size:	7850 byte(s)
Log Message:	MERGE OpenMD development 1783:1878 into trunk
#	Content
1	/*
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	}