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
#	User	Rev	Content
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			}