atomsGL.cpp 19.5 KB
Newer Older
1
2
3
4
5
#include <math.h>

#include "eprintf.h"
#include "TessShaders.h"
#include "UnitCellShaders.h"
6
#include "markerShaders.h"
7
8
9
10
11
#include "atomsGL.h"
#include "atoms.hpp"
#include "ConfigFile.h"
#include "CompileGLShader.h"
#include "polyhedron.h"
12
#include "Grid.h"
13

14
15
16
17
18
19
20
21
int getAtomTimesteps() 
{
	if (fixedAtoms)
		return 1;
	else
		return TIMESTEPS;
}

22
23
24
25
GLenum atomTexture(GLuint t)
{
	GLenum e;
	//rgh: scale atoms here
26
27
28
29
30
31
32
33
	//in google cardboard, this is called again if the program is running, so leave original or atoms get progresivelly smaller!
	float *a=new float[118*4];
	for (int i = 0; i < 118; i++) {
		a[i*4+0]=atomColours[i][0];
		a[i*4+1]=atomColours[i][1];
		a[i*4+2]=atomColours[i][2];
		a[i*4+3]=atomColours[i][3] * atomScaling;
	}
34
35
36
37
38
	glBindTexture(GL_TEXTURE_2D, t); //atom texture
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
39
	glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA32F, 118, 1, 0, GL_RGBA, GL_FLOAT, a);
40
41
42
43
44

	glBindTexture( GL_TEXTURE_2D, 0 );
	if ((e = glGetError()) != GL_NO_ERROR) {
		eprintf( "opengl error %d, atomTexture\n", e);
	}
45
	delete [] a;
46
47
48
49
50
51
	return e;
}

//WARNING: This should be called after SetupAtoms
//This means that numAtoms now has the cummulative distribution!
//This should be called after the atom texture is prepared, and therefore has the atomscaling pre-multiplied
52
53
54
GLenum SetupAtomsNoTess (GLuint **AtomVAO /*[3]*/, GLuint **AtomVertBuffer/*[2]*/, GLuint **AtomIndexBuffer/*[3]*/)
	//atoms, cloned atoms
	//rgh: FIXME: add AtomVAO[2] for atom trajectories
55
56
57
58
59
{
	//eprintf ("SetupAtomsNoTess 1");
if (!numAtoms)
		return 0;

60
61
62
63
64
if (!solid) {
	eprintf ("SetupAtomsNoTess, error: no solid defined");
	return 0;
}

65
66
67
68
69
70
//eprintf ("SetupAtomsNoTess 2");
	//for now, render an icosahedron
	//http://prideout.net/blog/?p=48 //public domain code
	//xyz nxnynz u=atom type ; 7 floats
	int e;

71
	int totalatoms=numAtoms[getAtomTimesteps() -1];
72
73
	
//eprintf ("SetupAtomsNoTess 2");
74
	*AtomVAO = new GLuint[3]; //atoms, cloned atoms, unused (bonds use Tess atom positions)
75
76
	*AtomIndexBuffer= new GLuint[3];//atoms, cloned atoms, bonds
	*AtomVertBuffer = new GLuint[2];//atoms, cloned atoms
77

78
	glGenVertexArrays(3, *AtomVAO);
79
80
81
82
83
84
85
86
87
88
89
	glGenBuffers(2, *AtomIndexBuffer);
	glGenBuffers(2, *AtomVertBuffer);
//eprintf ("SetupAtomsNoTess 3");
	glBindVertexArray((*AtomVAO)[0]);
	glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, (*AtomIndexBuffer)[0]);
	glBindBuffer(GL_ARRAY_BUFFER, (*AtomVertBuffer)[0]);
//eprintf ("SetupAtomsNoTess 4");
	glEnableVertexAttribArray(0);
	glEnableVertexAttribArray(1);
	glEnableVertexAttribArray(2);
	glDisableVertexAttribArray(3);
90
91
	//eprintf ("SetupAtomsNoTess 5, totalatoms=%d, nVerts=%d", totalatoms, solid->nVerts);
	float *tmp = new float[solid->nVerts * 7 * totalatoms];
92
93
	//eprintf ("SetupAtomsNoTess 6");
#ifdef INDICESGL32		
94
	int *tmpi = new int[solid->nFaces*3 * totalatoms];
95
96
97
	//eprintf ("SetupAtomsNoTess 7");
	int *currenti=tmpi;
#else
98
	unsigned short *tmpi = new unsigned short[solid->nFaces*3 * totalatoms];
99
100
101
102
103
104
	//eprintf ("SetupAtomsNoTess 7B");
	unsigned short *currenti=tmpi;
#endif

	float *current=tmp;
	//eprintf ("Before For 1");
105
	for (int p=0;p<getAtomTimesteps() ;p++) {
106
107
		for (int a = 0; a < numAtoms[p]-(p==0?0:numAtoms[p-1]); a++) {
			const int atomNumber = static_cast<int>(atoms[p][4 * a + 3]);
108
			const float radius = atomRadius(atomNumber)*atomScaling;
109
			for (int i = 0; i < solid->nVerts; i++) { //verts
110
				for (int k = 0; k < 3; k++) {
111
					*current++ = solid->Verts[3 * i + k]* radius +atoms[p][4 * a + k]; //pos
112
113
				}
				for (int k = 0; k < 3; k++) {
114
					*current++ = solid->Verts[3 * i + k]; //normal
115
				}
116
				*current++ = static_cast<float>(atomNumber);
117
			} //i
118
119
			for (int i = 0; i < solid->nFaces * 3; i++)
				*currenti++ = solid->Faces[i] + (a+(p==0?0:numAtoms[p-1]))*solid->nVerts;
120
121
		} //a
	} //p
122
		glBufferData(GL_ARRAY_BUFFER, sizeof(float) *totalatoms* 7 * solid->nVerts, tmp,
123
124
125
126
127
128
129
130
131
132
			GL_STATIC_DRAW);
		if ((e = glGetError()) != GL_NO_ERROR)
			eprintf("opengl error %d, glBufferData, l %d\n", e, __LINE__);

		glBufferData(GL_ELEMENT_ARRAY_BUFFER, 
#ifdef INDICESGL32		
	sizeof(int)
#else
	sizeof(unsigned int)
#endif
133
		* totalatoms * 3 * solid->nFaces, tmpi, GL_STATIC_DRAW);
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150


	glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 7*sizeof(float), (const void *)0);
	glVertexAttribPointer(1, 3, GL_FLOAT, GL_FALSE, 7 * sizeof(float), (const void *)(3*sizeof(float)));
	glVertexAttribPointer(2, 1, GL_FLOAT, GL_FALSE, 7 * sizeof(float), (const void *)(6 * sizeof(float)));

	if (glGetError() != GL_NO_ERROR)
		eprintf("opengl error attrib pointer 0\n");

	//glBindVertexArray(0);
	//glDisableVertexAttribArray(0);
	delete[] tmp;
	delete[] tmpi;
	if ((e = glGetError()) != GL_NO_ERROR)
		eprintf("opengl error %d, end of SetupAtoms, l %d\n", e, __LINE__);

	//FIXME TODO: cloned atoms
151
	tmp = new float[solid->nVerts * 7 * numClonedAtoms];
152
153
154
	current=tmp;
	//eprintf ("SetupAtomsNoTess 6");
#ifdef INDICESGL32		
155
	tmpi = new int[solid->nFaces*3 * numClonedAtoms];
156
157
158
	//eprintf ("SetupAtomsNoTess 7");
	currenti=tmpi;
#else
159
	tmpi = new unsigned short[solid->nFaces*3 * numClonedAtoms];
160
161
162
163
164
165
166
	//eprintf ("SetupAtomsNoTess 7B");
	currenti=tmpi;
#endif
	//eprintf ("Before For 2");

	for (int a = 0; a < numClonedAtoms; a++) {
		const int atomNumber = static_cast<int>(clonedAtoms[0][4 * a + 3]);
167
		const float radius = atomRadius(atomNumber)*atomScaling;
168
		for (int i = 0; i < solid->nVerts; i++) { //verts
169
				for (int k = 0; k < 3; k++) {
170
					*current++ = solid->Verts[3 * i + k]* radius +clonedAtoms[0][4 * a + k]; //pos
171
172
				}
				for (int k = 0; k < 3; k++) {
173
					*current++ = solid->Verts[3 * i + k]; //normal
174
				}
175
				*current++ =  static_cast<float>(atomNumber);
176
		} //i
177
178
		for (int i = 0; i < solid->nFaces * 3; i++)
			*currenti++ = solid->Faces[i] + a*solid->nVerts;
179
180
181
182
183
184
185
186
187
188
189
190
191
192
	} //a
	
	//eprintf ("After For 2");


	glBindVertexArray((*AtomVAO)[1]);
	glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, (*AtomIndexBuffer)[1]);
	glBindBuffer(GL_ARRAY_BUFFER, (*AtomVertBuffer)[1]);

	glEnableVertexAttribArray(0);
	glEnableVertexAttribArray(1);
	glEnableVertexAttribArray(2);
	glDisableVertexAttribArray(3);

193
	glBufferData(GL_ARRAY_BUFFER, sizeof(float) *numClonedAtoms* 7 * solid->nVerts, tmp,
194
195
196
197
198
199
200
201
202
203
			GL_STATIC_DRAW);
		if ((e = glGetError()) != GL_NO_ERROR)
			eprintf("opengl error %d, glBufferData, l %d\n", e, __LINE__);
	//eprintf ("After bufferdata, array buffer");
		glBufferData(GL_ELEMENT_ARRAY_BUFFER, 
#ifdef INDICESGL32		
	sizeof(int)
#else
	sizeof(unsigned int)
#endif
204
		* numClonedAtoms * 3 * solid->nFaces, tmpi, GL_STATIC_DRAW);
205
206
207
208
209
210
211
212
213
214
	//eprintf ("After bufferdata, element array buffer");
	glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 7*sizeof(float), (const void *)0);
	glVertexAttribPointer(1, 3, GL_FLOAT, GL_FALSE, 7 * sizeof(float), (const void *)(3*sizeof(float)));
	glVertexAttribPointer(2, 1, GL_FLOAT, GL_FALSE, 7 * sizeof(float), (const void *)(6 * sizeof(float)));

	if ((e = glGetError()) != GL_NO_ERROR)
		eprintf("opengl error %d, glVertexAttribPointer, l %d\n", e, __LINE__);

	delete[] tmp;
	delete[] tmpi;
215

216
	glBindVertexArray(0);
217
218
219
220
	return e;
} //SetupAtomsNoTess


221
GLenum SetupAtoms(GLuint **AtomVAO /*[3]*/, GLuint **AtomVertBuffer /*[2]*/, GLuint *BondIndices)
222
223
{
	if (!numAtoms)
224
		return glGetError();
225
226
227
228
229
230
231
	//rgh FIXME: put this all in the same vao
	
	//http://prideout.net/blog/?p=48 //public domain code
	//xyz u=atom type ; 4 floats
	int e;

	int totalatoms=0;
232
	for (int i=0;i<getAtomTimesteps() ;i++) {
233
234
235
236
		totalatoms += numAtoms[i];
	}
	eprintf("SetupAtoms: totalatoms=%d", totalatoms);

237
	*AtomVAO = new GLuint[3]; //atoms, cloned atoms, bonds //rgh fixme: for trajectories, we want to create another vao 
238
	*AtomVertBuffer = new GLuint[2];
239
240

	glGenVertexArrays(3, *AtomVAO);
241
	glGenBuffers(2, *AtomVertBuffer);
242
	glGenBuffers(1, BondIndices);
243
244
245
246
247
248
249
250

	glBindVertexArray((*AtomVAO)[0]);
	glBindBuffer(GL_ARRAY_BUFFER, (*AtomVertBuffer)[0]);

	glEnableVertexAttribArray(0);
	glEnableVertexAttribArray(1);
	glDisableVertexAttribArray(2);
	glDisableVertexAttribArray(3);
251
252
253
254

	e=glGetError();
	if (e!=GL_NO_ERROR)
		eprintf ("gl error %d, %s %d", e, __FILE__, __LINE__);
255
256
	float *tmp = new float[4 * totalatoms];
	float *current=tmp;
257
258
	
	const int atomlimit=30;
259
	const float bondscaling=0.7f;
260

261
262
	numBonds=new int[getAtomTimesteps() ];
	for (int p=0;p<getAtomTimesteps() ;p++) {
263

264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
			for (int a = 0; a < numAtoms[p]; a++) {
				for (int k = 0; k < 4; k++) {
					*current++ = atoms[p][4 * a + k];
				}
			} //a

		if (numAtoms[0]<atomlimit) {
		//eprintf ("searching bonds basic");
		//bonds FIXME quadractic complexity	
				for (int a1=0; a1 < numAtoms[p]; a1++) {
					for (int a2=a1+1; a2 < numAtoms[p]; a2++) {
						float d=0, r;
						for (int k=0;k<3;k++) {
							float dif=atoms[p][4 * a1 + k]-atoms[p][4 * a2 + k];
							d+=dif*dif;
						}
						r=atomRadius(static_cast<int>(atoms[p][4 * a1 + 3]))+
							atomRadius(static_cast<int>(atoms[p][4 * a2 + 3]));
282
						if (d*bondscaling<r*r) {// bond
283
284
285
286
287
288
289
290
291
292
293
							bonds.push_back(a1+(p==0?0:numAtoms[p-1]));
							bonds.push_back(a2+(p==0?0:numAtoms[p-1]));
						}
					}
				}
		} else { //more than 30 atoms, try grid optimization
		//eprintf ("searching bonds grid");

			float m[3];
			float M[3];
			for (int k=0; k<3;k++) {
294
				m[k]=M[k]=atoms[p][k];
295
			}
296
297
			for (int a = 1; a < numAtoms[p]; a++) {
				for (int k=0; k<3;k++) {
298
299
300
301
					if (m[k]>atoms[p][4*a+k])
						m[k]=atoms[p][4*a+k];
					if (M[k]<atoms[p][4*a+k])
						M[k]=atoms[p][4*a+k];
302
303
				}
			}
304
			grid g(m, M, pow(numAtoms[p], 1.0/3), bondscaling);
305
			for (int a = 1; a < numAtoms[p]; a++) 
306
				g.add(atoms[p]+4*a);
307
			for (int a = 0; a < numAtoms[p]; a++) {
308
				std::vector<float*> found=g.find(atoms[p]+4*a);
309
310
311
312
				for (int b=0;b<found.size();b++) {
					//if (found[b] < tmp+4*a) // already got this bound
					//	continue;
					bonds.push_back(a+(p==0?0:numAtoms[p-1]));
313
					bonds.push_back(((found[b]-atoms[p])/4)+(p==0?0:numAtoms[p-1]));
314
315
316
317
				}
			}
		}
		numBonds[p]=bonds.size();
318
319
320
		if (p!=0)
			numAtoms[p]+=numAtoms[p-1];
	} //p
321

322
323
324
	glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 4 * sizeof(float), (const void *)(0));
	glVertexAttribPointer(1, 1, GL_FLOAT, GL_FALSE, 4 * sizeof(float), (const void *)(3 * sizeof(float)));
	glBufferData(GL_ARRAY_BUFFER, sizeof(float) * totalatoms * 4 , tmp,
325
326
327
		GL_STATIC_DRAW);
	if ((e = glGetError()) != GL_NO_ERROR)
		eprintf( "opengl error %d, glBufferData, l %d\n", e, __LINE__);
328

329
	glBindVertexArray(0);
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347

	if ((e = glGetError()) != GL_NO_ERROR)
		eprintf( "opengl error %d, end of SetupAtoms, l %d\n", e, __LINE__);

	if (showTrajectories) {
			//fill the restart buffer
		//use abc for measuring
		float max=0;
		if (has_abc) {
			for (int i=0;i<3;i++)
				for (int j=0;j<3;j++)
				max+=abc[i][j];
			max /=9*2;
		}

		for (unsigned int t=0;t<atomtrajectories.size();t++) {
			atomtrajectoryrestarts.push_back(std::vector<int>());
			atomtrajectoryrestarts[t].push_back(0);
348
			for (int p=1;p<getAtomTimesteps() ;p++) {
349
				int a=atomtrajectories[t];
350
351
352
353
354
				if (has_abc)
					if (fabs(atoms[p][a*4+0]-atoms[p-1][a*4+0])+
						fabs(atoms[p][a*4+1]-atoms[p-1][a*4+1])+
						fabs(atoms[p][a*4+2]-atoms[p-1][a*4+2])>max)
							atomtrajectoryrestarts[t].push_back(p);
355
			}
356
			atomtrajectoryrestarts[t].push_back(getAtomTimesteps() );
357
358
359
		}
	}
	delete[] tmp;
360
361
362
363
364
365
366
367
368
369
370
371
	//bonds
	glBindVertexArray((*AtomVAO)[2]);
	glBindBuffer(GL_ARRAY_BUFFER, (*AtomVertBuffer)[0]);
	glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, *BondIndices);
	glBufferData(GL_ELEMENT_ARRAY_BUFFER, sizeof(int)*bonds.size(), bonds.data(), GL_STATIC_DRAW);
	glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 4 * sizeof(float), (const void *)(0));
	glEnableVertexAttribArray(0);
	glBindVertexArray(0);

	e=glGetError();
	if ((e = glGetError()) != GL_NO_ERROR)
		eprintf( "opengl error %d, creating chemical bonds, l %d\n", e, __LINE__);
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391

	//now clones
	if (basisvectorreps ||!clonedAtoms) //do not replicate
		return e;


	glBindVertexArray((*AtomVAO)[1]); //rgh FIXME, only works for TIMESTEPS=1
	glBindBuffer(GL_ARRAY_BUFFER, (*AtomVertBuffer)[1]);
	glBufferData(GL_ARRAY_BUFFER, sizeof(float) * clonedAtoms[0].size(), clonedAtoms[0].data(),
			GL_STATIC_DRAW);
	glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 4 * sizeof(float), (const void *)(0));
	glVertexAttribPointer(1, 1, GL_FLOAT, GL_FALSE, 4 * sizeof(float), (const void *)(3 * sizeof(float)));
	glEnableVertexAttribArray(0);
	glEnableVertexAttribArray(1);
	if ((e = glGetError()) != GL_NO_ERROR)
		eprintf( "opengl error %d, end of Setup cloned Atoms, l %d\n", e, __LINE__);

	//rgh: we will need these again if we don't have tesselation
	//delete[] clonedAtoms;
	//clonedAtoms=0;
392
	glBindVertexArray(0);
393
394
395
	return e;
}

396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
GLenum SetupMarker(GLuint *MarkerVAO, GLuint *MarkerVertBuffer)
{
	if (!markers)
		return glGetError();
	GLenum e;
	if ((e = glGetError()) != GL_NO_ERROR)
		eprintf( "opengl error %d, begin of SetupMarker\n", e, __LINE__);

	glGenVertexArrays(1, MarkerVAO);
	glGenBuffers(1, MarkerVertBuffer);
	if ((e = glGetError()) != GL_NO_ERROR)
		eprintf( "opengl error %d, glGenBuffers, l %d\n", e, __LINE__);

	glBindVertexArray(*MarkerVAO);
	glBindBuffer(GL_ARRAY_BUFFER, *MarkerVertBuffer);

	glEnableVertexAttribArray(0);
	glEnableVertexAttribArray(1);
	glDisableVertexAttribArray(2);
	glDisableVertexAttribArray(3);

	const float size=atomRadius(0)*atomScaling*markerscaling;
	float *tmp = new float [8*TIMESTEPS];
	for (int i=0;i<TIMESTEPS;i++) {
		for (int j=0;j<3;j++) { //center [3]
			tmp[i*8+j]=markers[i][j];
		}
		tmp[i*8+3]=0.8*size; //size [1]
		for (int j=0;j<4;j++) {//colour[4]
			tmp[i*8+4+j]=markercolours[i][j];
		}
	}
	
	glBufferData(GL_ARRAY_BUFFER, sizeof(float) * TIMESTEPS*8 , tmp,
			GL_STATIC_DRAW);
	glVertexAttribPointer(0, 4, GL_FLOAT, GL_FALSE, 8 * sizeof(float), (const void *)(0));
	glVertexAttribPointer(1, 4, GL_FLOAT, GL_FALSE, 8 * sizeof(float), (const void *)(4*sizeof(float)));
	glBindVertexArray(0);
	return glGetError();
}

437
438
GLenum SetupUnitCell(GLuint *UnitCellVAO, GLuint *UnitCellVertBuffer, GLuint *UnitCellIndexBuffer)
{
439
	//add here both unit cell and supercell
440
	GLenum e;
441
442
443
444
	if ((e = glGetError()) != GL_NO_ERROR)
		eprintf( "opengl error %d, begin of SetupUnitCell\n", e, __LINE__);
	if (!has_abc)
		return e;
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
	glGenVertexArrays(1, UnitCellVAO);
	glGenBuffers(1, UnitCellVertBuffer);
	glGenBuffers(1, UnitCellIndexBuffer);
	if ((e = glGetError()) != GL_NO_ERROR)
		eprintf( "opengl error %d, glGenBuffers, l %d\n", e, __LINE__);

	glBindVertexArray(*UnitCellVAO);
	glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, *UnitCellIndexBuffer);
	glBindBuffer(GL_ARRAY_BUFFER, *UnitCellVertBuffer);

	glEnableVertexAttribArray(0);
	glDisableVertexAttribArray(1);
	glDisableVertexAttribArray(2);
	glDisableVertexAttribArray(3);

460
	float *tmp = new float[3*8*2];
461
	//0, a, b, c, a+b+c, b+c, a+c, a+b
462
	for (int i=0;i<3;i++) { //unit cell
463
464
465
466
467
468
469
470
		tmp[0+i]=0;
		for (int j=0;j<3;j++)
			tmp[3*(j+1)+i]=abc[j][i];
		tmp[3*4+i]=abc[0][i]+abc[1][i]+abc[2][i];
		tmp[3*5+i]=			abc[1][i]+abc[2][i];
		tmp[3*6+i]=abc[0][i]+		abc[2][i];
		tmp[3*7+i]=abc[0][i]+abc[1][i];
	}
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
	float displ[3]={0,0,0};
	if (translations && ISOS) 
		for (int j=0;j<3;j++)
			for (int i=0;i<3;i++)
				displ[i]+=-translations[0][i]*abc[j][i];
	for (int i=0;i<3;i++) { //rgh fixme, add displacement here as well
		tmp[3*8+i]=displ[i];
		for (int j=0;j<3;j++)
			tmp[3*(j+8+1)+i]=abc[j][i]*supercell[j]+displ[i];
		tmp[3*12+i]=abc[0][i]*supercell[0]+abc[1][i]*supercell[1]+abc[2][i]*supercell[2]+displ[i];
		tmp[3*13+i]=			abc[1][i]*supercell[1]+abc[2][i]*supercell[2]+displ[i];
		tmp[3*14+i]=abc[0][i]*supercell[0]+		abc[2][i]*supercell[2]+displ[i];
		tmp[3*15+i]=abc[0][i]*supercell[0]+abc[1][i]*supercell[1]+displ[i];
	}
	int tmpi[12*2*2]={ //lines, unit cell, 
486
487
488
489
490
491
492
493
494
495
496
		0,1, 
		1,6,
		6,3,
		3,0,
		2,7,
		7,4,
		4,5,
		5,2,
		0,2,
		1,7,
		6,4,
497
498
499
500
501
502
503
504
505
506
507
508
509
		3,5, // supercell
		0+8,1+8, 
		1+8,6+8,
		6+8,3+8,
		3+8,0+8,
		2+8,7+8,
		7+8,4+8,
		4+8,5+8,
		5+8,2+8,
		0+8,2+8,
		1+8,7+8,
		6+8,4+8,
		3+8,5+8
510
	};
511
	glBufferData(GL_ARRAY_BUFFER, sizeof(float) * 3*8*2 , tmp,
512
513
			GL_STATIC_DRAW);
	if ((e = glGetError()) != GL_NO_ERROR)
514
		eprintf( "opengl error %d, glBufferData vertex, l %d\n", e, __LINE__);
515
516
	glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 3 * sizeof(float), (const void *)(0));
	glBufferData(GL_ELEMENT_ARRAY_BUFFER, sizeof(tmpi), tmpi, GL_STATIC_DRAW);
517
518
	if ((e = glGetError()) != GL_NO_ERROR)
		eprintf( "opengl error %d, glBufferData index, l %d\n", e, __LINE__);
519
	glBindVertexArray(0);
520
521
522
523
	return e;
}


524
525
526
bool PrepareUnitCellAtomShader (GLuint *AtomP, GLuint *cellP, GLuint *MarkerP, 
								GLint *AtomMatrixLocation, GLint *UnitCellMatrixLocation,  GLint *UnitCellColourLocation,
								GLint *MarkerMatrixLocation){
527
528
529
530
531
532
	if (!PrepareAtomShader(AtomP, AtomMatrixLocation))
		return false;

	if (!PrepareUnitCellShader(cellP, UnitCellMatrixLocation, UnitCellColourLocation))
		return false;

533
534
535
	if (!PrepareMarkerShader(MarkerP, MarkerMatrixLocation))
		return false;

536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
	return true;
}

bool PrepareAtomShader (GLuint *AtomP, GLint *AtomMatrixLocation){
		//https://www.gamedev.net/topic/591110-geometry-shader-point-sprites-to-spheres/
	//no rotation, only translations means we can do directional lighting in the shader.
	//FIXME
	//http://stackoverflow.com/questions/40101023/flat-shading-in-webgl
	*AtomP = CompileGLShader(
		AtomShaders[SHADERNAME],
		AtomShaders[SHADERVERTEX],
		AtomShaders[SHADERFRAGMENT],
		AtomShaders[SHADERTESSEVAL]
		);
	*AtomMatrixLocation=glGetUniformLocation(*AtomP, "matrix");
	if( *AtomMatrixLocation == -1 )
	{
		eprintf( "Unable to find matrix uniform in atom shader\n" );
		return false;
	}
	return true;
}

559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
bool PrepareMarkerShader (GLuint *MP, GLint *MMatrixLocation){
		//https://www.gamedev.net/topic/591110-geometry-shader-point-sprites-to-spheres/
	//no rotation, only translations means we can do directional lighting in the shader.
	//FIXME
	//http://stackoverflow.com/questions/40101023/flat-shading-in-webgl
	*MP = CompileGLShader(
		MarkerShaders[SHADERNAME],
		MarkerShaders[SHADERVERTEX],
		MarkerShaders[SHADERFRAGMENT],
		MarkerShaders[SHADERTESSEVAL]
		);
	*MMatrixLocation=glGetUniformLocation(*MP, "matrix");
	if( *MMatrixLocation == -1 )
	{
		eprintf( "Unable to find matrix uniform in atom shader\n" );
		return false;
	}
	return true;
}

579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
bool PrepareAtomShaderNoTess (GLuint *AtomP, GLint *AtomMatrixLocation){
		//https://www.gamedev.net/topic/591110-geometry-shader-point-sprites-to-spheres/
	//no rotation, only translations means we can do directional lighting in the shader.
	//FIXME
	//http://stackoverflow.com/questions/40101023/flat-shading-in-webgl
	*AtomP = CompileGLShader(
		AtomShadersNoTess[SHADERNAME],
		AtomShadersNoTess[SHADERVERTEX],
		AtomShadersNoTess[SHADERFRAGMENT],
		AtomShadersNoTess[SHADERTESSEVAL]
		);
	*AtomMatrixLocation=glGetUniformLocation(*AtomP, "matrix");
	if( *AtomMatrixLocation == -1 )
	{
		eprintf( "Unable to find matrix uniform in atom shader no tess\n" );
		return false;
	}
	return true;
}


bool PrepareUnitCellShader (GLuint *cellP, GLint *UnitCellMatrixLocation,  GLint *UnitCellColourLocation){
	*cellP= CompileGLShader(
		UnitCellShaders[SHADERNAME],
		UnitCellShaders[SHADERVERTEX],
		UnitCellShaders[SHADERFRAGMENT],
		UnitCellShaders[SHADERTESSEVAL]
		);
	*UnitCellMatrixLocation=glGetUniformLocation(*cellP, "matrix");
	if( *UnitCellMatrixLocation == -1 )
	{
		eprintf( "Unable to find matrix uniform in UnitCell shader\n" );
		return false;
	}
	*UnitCellColourLocation=glGetUniformLocation(*cellP, "color");
	if( *UnitCellColourLocation == -1 )
	{
		eprintf( "Unable to find color uniform in UnitCell shader\n" );
		return false;
	}
	return true;
}


/**p: input, f: output*/
void GetDisplacement(int p[3], float f[3])
{
float delta[3][3];
for (int ss=0;ss<3;ss++)
	for (int i=0;i<3;i++)
		delta[ss][i]=static_cast<float>(p[ss])*abc[ss][i];

for (int i=0;i<3;i++)
	f[i]=0;

for (int ss=0;ss<3;ss++)
	for (int i=0;i<3;i++)
		f[i]+=delta[ss][i];
637
}