1 /* Ergo, version 3.8, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2019 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4 * and Anastasia Kruchinina.
5 *
6 * This program is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 *
19 * Primary academic reference:
20 * Ergo: An open-source program for linear-scaling electronic structure
21 * calculations,
22 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23 * Kruchinina,
24 * SoftwareX 7, 107 (2018),
25 * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26 *
27 * For further information about Ergo, see <http://www.ergoscf.org>.
28 */
29
30 #include <cstdio>
31 #include <cmath>
32 #include <cstdlib>
33 #include <string>
34 #include <cstring>
35 #include <vector>
36 #include <cassert>
37
38 typedef struct {
39 std::string label;
40 double coord[3];
41 } atom_struct;
42
43 struct vec3d_struct {
44 double coord[3];
setvec3d_struct45 void set(double x, double y, double z) {
46 coord[0] = x;
47 coord[1] = y;
48 coord[2] = z;
49 }
vec3d_structvec3d_struct50 vec3d_struct(double x, double y, double z) {
51 set(x, y, z);
52 }
53 };
54
write_xyz_file(int nAtoms,atom_struct atomList[],const char * outFileName,const char * comment)55 int write_xyz_file(int nAtoms, atom_struct atomList[], const char* outFileName, const char* comment) {
56 const int maxLabelLen = 4;
57 char paddingString[888];
58 int k, len;
59 FILE* f = fopen(outFileName, "wt");
60 if(f == NULL) {
61 printf("error opening file '%s' for writing\n", outFileName);
62 return -1;
63 }
64 // First line should cantain the number of atoms
65 fprintf(f, "%i\n", nAtoms);
66 // Second line is a comment line
67 fprintf(f, "%s\n", comment);
68 for(int i = 0; i < nAtoms; i++) {
69 len = strlen(atomList[i].label.c_str());
70 if(len > maxLabelLen) {
71 printf("error: label too long: '%s'\n", atomList[i].label.c_str());
72 return -1;
73 }
74 for(k = 0; k < maxLabelLen-len; k++)
75 paddingString[k] = ' ';
76 paddingString[k] = '\0';
77 fprintf(f, "%s%s %15.9f %15.9f %15.9f\n",
78 atomList[i].label.c_str(),
79 paddingString,
80 atomList[i].coord[0],
81 atomList[i].coord[1],
82 atomList[i].coord[2]);
83 } // END FOR i
84 fclose(f);
85 printf("File '%s' created OK\n", outFileName);
86 return 0;
87 }
88
get_dist(vec3d_struct const & pos1,vec3d_struct const & pos2)89 static double get_dist(vec3d_struct const & pos1, vec3d_struct const & pos2) {
90 double sum_of_squares = 0;
91 for(int i = 0; i < 3; i++) {
92 double d = pos1.coord[i] - pos2.coord[i];
93 sum_of_squares += d*d;
94 }
95 return sqrt(sum_of_squares);
96 }
97
get_vec_length(vec3d_struct const & A)98 static double get_vec_length(vec3d_struct const & A) {
99 double sum = 0;
100 for(int i = 0; i < 3; i++)
101 sum += A.coord[i]*A.coord[i];
102 return sqrt(sum);
103 }
104
get_dot_product(vec3d_struct const & A,vec3d_struct const & B)105 static double get_dot_product(vec3d_struct const & A, vec3d_struct const & B) {
106 double sum = 0;
107 for(int i = 0; i < 3; i++)
108 sum += A.coord[i]*B.coord[i];
109 return sum;
110 }
111
normalize_vector(vec3d_struct & A)112 static void normalize_vector(vec3d_struct & A) {
113 double d = get_vec_length(A);
114 for(int i = 0; i < 3; i++)
115 A.coord[i] /= d;
116 }
117
add_to_vec(vec3d_struct & A,vec3d_struct const & B,double factor)118 static void add_to_vec(vec3d_struct & A, vec3d_struct const & B, double factor) {
119 for(int i = 0; i < 3; i++)
120 A.coord[i] += factor * B.coord[i];
121 }
122
main(int argc,char * argv[])123 int main(int argc, char* argv[])
124 {
125 printf("generate_alkane version 1.0\n");
126 if(argc != 3) {
127 printf("Usage:\n");
128 printf("generate_alkane N filename\n");
129 printf("where\n");
130 printf(" N : number of C atoms.\n");
131 printf(" filename : filename of new xyz file to create.\n");
132 return -1;
133 }
134 int N = atoi(argv[1]);
135 if(N < 1) {
136 printf("Error: (N < 1)\n");
137 return -1;
138 }
139 const char* filename = argv[2];
140 printf("N = %d\n", N);
141 printf("filename = '%s'\n", filename);
142
143 const double C_C_dist = 1.54;
144 const double C_H_dist = 1.09;
145
146 // To get the angles right, start by setting up unit vectors pointing in the directions we need.
147 // Create perfect tetrahedron A-B-C-D
148 double h = 1 / sqrt(8);
149 vec3d_struct A(0, 1, -h);
150 vec3d_struct B(sqrt(3)/2, -0.5, -h);
151 vec3d_struct C(-sqrt(3)/2, -0.5, -h);
152 vec3d_struct D(0, 0, sqrt(1+h*h));
153 // Check that tetrahedron ABCD is correct; distance to origin should be same for all of the points, and the distances between them should also be the same.
154 double lengthA = get_vec_length(A);
155 double lengthB = get_vec_length(B);
156 double lengthC = get_vec_length(C);
157 double lengthD = get_vec_length(D);
158 double tol = 1e-8;
159 assert(fabs(lengthA-lengthB) < tol);
160 assert(fabs(lengthA-lengthC) < tol);
161 assert(fabs(lengthA-lengthD) < tol);
162 double distAB = get_dist(A, B);
163 double distAC = get_dist(A, C);
164 double distAD = get_dist(A, D);
165 double distBC = get_dist(B, C);
166 double distBD = get_dist(B, D);
167 double distCD = get_dist(C, D);
168 assert(fabs(distAB-distAC) < tol);
169 assert(fabs(distAB-distAD) < tol);
170 assert(fabs(distAB-distBC) < tol);
171 assert(fabs(distAB-distBD) < tol);
172 assert(fabs(distAB-distCD) < tol);
173 // Normalize vectors
174 normalize_vector(A);
175 normalize_vector(B);
176 normalize_vector(C);
177 normalize_vector(D);
178
179 // Check angle between vectors
180 double dot_product_AB = get_dot_product(A, B);
181 double dot_product_AC = get_dot_product(A, C);
182 double dot_product_AD = get_dot_product(A, D);
183 double dot_product_BC = get_dot_product(B, C);
184 double dot_product_BD = get_dot_product(B, D);
185 double dot_product_CD = get_dot_product(C, D);
186 assert(fabs(dot_product_AB-dot_product_AC) < tol);
187 assert(fabs(dot_product_AB-dot_product_AD) < tol);
188 assert(fabs(dot_product_AB-dot_product_BC) < tol);
189 assert(fabs(dot_product_AB-dot_product_BD) < tol);
190 assert(fabs(dot_product_AB-dot_product_CD) < tol);
191 double angle = acos(dot_product_AB);
192 double pi = 2 * asin(1);
193 double angle_in_degrees = angle * 180 / pi;
194 printf("Angle between vectors: %15.10f <--> %15.10f degrees\n", angle, angle_in_degrees);
195
196 int NC = N;
197
198 // Place C atoms
199 std::vector<vec3d_struct> C_pts;
200 C_pts.reserve(NC);
201
202 vec3d_struct currPos(0, 0, 0);
203 for(int i = 0; i < NC; i++) {
204 C_pts.push_back(currPos);
205 if(i % 2 == 0) {
206 // Move in direction D
207 add_to_vec(currPos, D, C_C_dist);
208 }
209 else {
210 // Move in direction -A
211 add_to_vec(currPos, A, -C_C_dist);
212 }
213 }
214
215 // Now add H atoms for each C atom
216 std::vector<vec3d_struct> H_pts;
217 H_pts.reserve(2*NC+2);
218 for(int i = 0; i < NC; i++) {
219 vec3d_struct Cpos = C_pts[i];
220 if(i == 0) {
221 // Special case: first C atom, should have 3 H atoms.
222 vec3d_struct Hpos1 = Cpos;
223 add_to_vec(Hpos1, A, C_H_dist);
224 H_pts.push_back(Hpos1);
225 vec3d_struct Hpos2 = Cpos;
226 add_to_vec(Hpos2, B, C_H_dist);
227 H_pts.push_back(Hpos2);
228 vec3d_struct Hpos3 = Cpos;
229 add_to_vec(Hpos3, C, C_H_dist);
230 H_pts.push_back(Hpos3);
231 if(NC == 1) {
232 // Add 4th H atom in this case
233 vec3d_struct Hpos4 = Cpos;
234 add_to_vec(Hpos4, D, C_H_dist);
235 H_pts.push_back(Hpos4);
236 }
237 }
238 else if(i == NC-1) {
239 // Special case: last C atom, should have 3 H atoms.
240 if(i % 2 == 0) {
241 vec3d_struct Hpos1 = Cpos;
242 add_to_vec(Hpos1, D, C_H_dist);
243 H_pts.push_back(Hpos1);
244 vec3d_struct Hpos2 = Cpos;
245 add_to_vec(Hpos2, B, C_H_dist);
246 H_pts.push_back(Hpos2);
247 vec3d_struct Hpos3 = Cpos;
248 add_to_vec(Hpos3, C, C_H_dist);
249 H_pts.push_back(Hpos3);
250 }
251 else {
252 vec3d_struct Hpos1 = Cpos;
253 add_to_vec(Hpos1, A, -C_H_dist);
254 H_pts.push_back(Hpos1);
255 vec3d_struct Hpos2 = Cpos;
256 add_to_vec(Hpos2, B, -C_H_dist);
257 H_pts.push_back(Hpos2);
258 vec3d_struct Hpos3 = Cpos;
259 add_to_vec(Hpos3, C, -C_H_dist);
260 H_pts.push_back(Hpos3);
261 }
262 }
263 else {
264 // Normal case, 2 H atoms.
265 if(i % 2 == 0) {
266 vec3d_struct Hpos1 = Cpos;
267 add_to_vec(Hpos1, B, C_H_dist);
268 H_pts.push_back(Hpos1);
269 vec3d_struct Hpos2 = Cpos;
270 add_to_vec(Hpos2, C, C_H_dist);
271 H_pts.push_back(Hpos2);
272 }
273 else {
274 vec3d_struct Hpos1 = Cpos;
275 add_to_vec(Hpos1, B, -C_H_dist);
276 H_pts.push_back(Hpos1);
277 vec3d_struct Hpos2 = Cpos;
278 add_to_vec(Hpos2, C, -C_H_dist);
279 H_pts.push_back(Hpos2);
280 }
281 }
282 }
283 int NH = H_pts.size();
284 assert(NH == 2*NC+2);
285
286 // Now create final list of atoms
287 int nAtomsTot = NC + NH;
288 std::vector<atom_struct> atoms;
289 atoms.reserve(nAtomsTot);
290 // Add C atoms
291 for(int i = 0; i < NC; i++) {
292 atom_struct newAtom;
293 newAtom.label = "C";
294 for(int k = 0; k < 3; k++)
295 newAtom.coord[k] = C_pts[i].coord[k];
296 atoms.push_back(newAtom);
297 }
298 // Add H atoms
299 for(int i = 0; i < NH; i++) {
300 atom_struct newAtom;
301 newAtom.label = "H";
302 for(int k = 0; k < 3; k++)
303 newAtom.coord[k] = H_pts[i].coord[k];
304 atoms.push_back(newAtom);
305 }
306
307 char comment[888];
308 sprintf(comment, "C%dH%d xyz file created by the generate_alkane program, NC = %d, NH = %d, C_C_dist = %f, C_H_dist = %f", NC, NH, NC, NH, C_C_dist, C_H_dist);
309 int nAtoms = atoms.size();
310 if(write_xyz_file(nAtoms, &atoms[0], filename, comment) != 0) {
311 printf("Error in write_xyz_file().\n");
312 return -1;
313 }
314
315 printf("Done. Alkane xyz file '%s' created, %d C atoms and %d H atoms, %d atoms in total.\n", filename, NC, NH, nAtoms);
316
317 return 0;
318 }
319