Allow numbered atoms and unknown atom types in xyz

Atomic bonds: optimized creation with 3D grid
Fix atom trajectories if no cell vectors available

NOMADVRLib/ConfigFile.cpp

@@ -21,6 +21,8 @@ int* numAtoms; //[timesteps]
 float **atoms; //[timesteps][numAtoms[i]*4] //xyzu, u=atom number
 float atomScaling;
 std::vector<float> *clonedAtoms;
+std::vector<int> bonds;
+int *numBonds;
 int numClonedAtoms;
 int *basisvectorreps;
 
@@ -127,6 +129,7 @@ int loadConfigFile(const char * f)
 	numAtoms=0;
 	atomScaling=1;
 	clonedAtoms=0;
+	bonds.clear();
 	showTrajectories = false;
 	basisvectorreps=0;
 	numClonedAtoms=0;
@@ -396,6 +399,11 @@ int loadConfigFile(const char * f)
 		for (int i=0;i<*numAtoms;i++)
 			atomtrajectories.push_back(i);
 	}
+
+	//chemical bonds
+	//if (numAtoms) {
+	//}
+
 	//eprintf ("Before returning, numatoms 0 =%d", numAtoms[0]);
 	return 0;
 }

NOMADVRLib/ConfigFile.h

@@ -17,6 +17,8 @@ extern int* numAtoms; //[timesteps]
 extern float **atoms; //[timesteps][numAtoms[i]*4] //xyzu, u=atom number
 extern float atomScaling;
 extern std::vector<float> *clonedAtoms;
+extern std::vector<int> bonds;
+extern int *numBonds;
 extern int numClonedAtoms;
 extern int *basisvectorreps;
 

NOMADVRLib/Grid.cpp

+#include "Grid.h"
+#include "atoms.hpp" //for radius
+
+grid::grid (float *m, float *M, int dims, float s):scale(s) {
+	content = new std::vector<float*> [dims*dims*dims];
+	for (int i=0;i<3;i++) {
+		this->m[i]=m[i];
+		this->M[i]=M[i];
+		dif[i]=M[i]-m[i];
+	}
+	this->dims=dims;
+	maxradius=0;
+}
+
+grid::~grid()
+{
+	delete [] content;
+}
+
+void grid::coordinates(const float pos[3], int c[3])
+{
+	for (int i=0;i<3;i++) {
+		c[i]=floor((pos[i]-m[i])/dif[i]*dims);
+		if (c[i]>=dims)
+			c[i]=dims-1;
+		else if (c[i]<0)
+			c[i]=0;
+	}
+}
+
+void grid::add (float *p) //compatible with the atoms xyzr
+{
+	int pos[3];
+	coordinates (p, pos);
+	
+	content[pos[0]*dims*dims + pos[1]*dims+pos[2]].push_back(p);
+	float ar=atomRadius(p[3]);
+	if (ar>maxradius)
+		maxradius=ar;
+}
+
+bool grid::compare (float *a, float *b)
+{
+	if (a<=b) //already returned when searching a beforehand
+		return false;
+	float d=0;
+	for (int i=0;i<3;i++)
+		d+=(a[i]-b[i])*(a[i]-b[i]);
+	if (d*scale < (a[3]+b[3]))
+		return true;
+}
+
+std::vector<float*> grid::find (float *p) 
+{
+	std::vector<float*> result;
+	//search a sphere centered in p, of radius (p[3]+maxradius) *scale
+
+	//start by searching a cube
+	int mc[3];
+	int Mc[3];
+	float mp[3];
+	float Mp[3];
+	for (int i=0;i<3;i++) {
+		mp[i]=p[i]-(atomRadius(p[3])+maxradius)/scale;
+		Mp[i]=p[i]+(atomRadius(p[3])+maxradius)/scale;
+	}
+	coordinates(mp, mc);
+	coordinates(Mp, Mc);
+	for (int x=mc[0];x<Mc[0];x++)
+		for (int y=mc[1];y<Mc[1];y++)
+			for (int z=mc[2];z<Mc[2];z++) {
+				const int c=x*dims*dims + y*dims+z;
+				for (int i=0;i<content[c].size();i++) {
+					if (compare (content[c][i], p))
+						result.push_back(content[c][i]);
+				}
+			}
+	return result;			
+}
+

NOMADVRLib/Grid.h

+#include <vector>
+
+class grid 
+{
+public:
+	grid (float *m, float *M, int dims, float s);
+	~grid();
+	void add (float *p); //compatible with the atoms xyzr
+	std::vector<float*> find (float *p);
+private:
+	void coordinates(const float pos[3], int c[3]);
+	bool compare (float *a, float *b);
+	std::vector <float*> *content;
+	float m[3];
+	float M[3];
+	float dif[3];
+	int dims;
+	float maxradius;
+	const float scale;
+}; // grid

NOMADVRLib/IsosurfacesGL.cpp

@@ -294,8 +294,11 @@ bool SetupDepthPeeling(int renderWidth, int renderHeight, int zlayers, GLuint *t
 			eprintf("opengl error %d, SetupDepthPeeling b\n", e);
 
 		//cleaned at each frame
+#if defined(WIN32) || defined(CAVE)
+		CleanDepthTexture(textures[i]);
+#else
 		CleanDepthTexture(textures[i], renderWidth, renderHeight);
-
+#endif
 		if ((e = glGetError()) != GL_NO_ERROR)
 			eprintf("opengl error %d, SetupDepthPeeling c\n", e);
 	}

NOMADVRLib/atoms.cpp

@@ -163,6 +163,11 @@ float atomColours[][4] =
 {MISSINGR, MISSINGG, MISSINGB, MISSINGRADIUS},//Og
 };
 
+float atomRadius (int i)
+{
+	return atomColours[i][3];
+}
+
 void discardline (FILE *F)
 {
 int c;
@@ -173,9 +178,19 @@ do {
 
 int findAtom(const char *const s)
 {
+	//discard number at end
+	char x[10];
+	char *p=x; 
+	const char *q=s;
+	while (*q!=0) {
+		if (*q >='0' && *q <= '9')
+			break;
+		*p++=*q++;
+	}
+	*p=0;
 	//rgh FIXME, add caching
 	for (unsigned int i=0;i<sizeof(atomNames)/sizeof(const char *);i++)
-		if (!strcmp(s, atomNames[i]))
+		if (!strcmp(x, atomNames[i]))
 			return i;
 	return -1;
 }
@@ -215,8 +230,11 @@ int readAtomsXYZ(const char *const file, int **numatoms, int *timesteps, float *
 			if (r<4)
 				return -2;
 			int a=findAtom(s);
-			if (a==-1)
-				return -3;
+			if (a==-1) {
+				eprintf ("Read atoms xyz, atom type unknown: %s", s);
+				a=83; //rgh fixme, this is unknown, unknown unknown, make a real unknown at the end.
+				//return -3;
+			}
 			(mypos.back())[4*i+3]=(float)a;
 		}
 	}

NOMADVRLib/atoms.hpp

@@ -29,6 +29,8 @@ extern const char * readAtomsXYZErrors[];
 extern const char * readAtomsCubeErrors[];
 extern const char * readAtomsJsonErrors[];
 
+float atomRadius (int i);
+
 //internal functions
 void discardline (FILE *F);
 void Clone (float tmppos[3], float k, std::vector<float>* clonedAtoms);

NOMADVRLib/atomsGL.cpp

@@ -8,6 +8,7 @@
 #include "ConfigFile.h"
 #include "CompileGLShader.h"
 #include "polyhedron.h"
+#include "Grid.h"
 
 GLenum atomTexture(GLuint t)
 {
@@ -39,7 +40,9 @@ GLenum atomTexture(GLuint t)
 //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
-GLenum SetupAtomsNoTess (GLuint **AtomVAO, GLuint **AtomVertBuffer, GLuint **AtomIndexBuffer)
+GLenum SetupAtomsNoTess (GLuint **AtomVAO /*[3]*/, GLuint **AtomVertBuffer/*[2]*/, GLuint **AtomIndexBuffer/*[3]*/)
+	//atoms, cloned atoms
+	//rgh: FIXME: add AtomVAO[2] for atom trajectories
 {
 	//eprintf ("SetupAtomsNoTess 1");
 if (!numAtoms)
@@ -59,11 +62,11 @@ if (!solid) {
 	int totalatoms=numAtoms[TIMESTEPS-1];
 	
 //eprintf ("SetupAtomsNoTess 2");
-	*AtomVAO = new GLuint[2]; //atoms, cloned atoms
-	*AtomIndexBuffer= new GLuint[2];
-	*AtomVertBuffer = new GLuint[2];
+	*AtomVAO = new GLuint[3]; //atoms, cloned atoms, bonds
+	*AtomIndexBuffer= new GLuint[3];//atoms, cloned atoms, bonds
+	*AtomVertBuffer = new GLuint[2];//atoms, cloned atoms
 
-	glGenVertexArrays(2, *AtomVAO);
+	glGenVertexArrays(3, *AtomVAO);
 	glGenBuffers(2, *AtomIndexBuffer);
 	glGenBuffers(2, *AtomVertBuffer);
 //eprintf ("SetupAtomsNoTess 3");
@@ -93,7 +96,7 @@ if (!solid) {
 	for (int p=0;p<TIMESTEPS;p++) {
 		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]);
-			const float radius = atomColours[atomNumber][3]*atomScaling;
+			const float radius = atomRadius(atomNumber)*atomScaling;
 			for (int i = 0; i < solid->nVerts; i++) { //verts
 				for (int k = 0; k < 3; k++) {
 					*current++ = solid->Verts[3 * i + k]* radius +atoms[p][4 * a + k]; //pos
@@ -152,7 +155,7 @@ if (!solid) {
 
 	for (int a = 0; a < numClonedAtoms; a++) {
 		const int atomNumber = static_cast<int>(clonedAtoms[0][4 * a + 3]);
-		const float radius = atomColours[atomNumber][3]*atomScaling;
+		const float radius = atomRadius(atomNumber)*atomScaling;
 		for (int i = 0; i < solid->nVerts; i++) { //verts
 				for (int k = 0; k < 3; k++) {
 					*current++ = solid->Verts[3 * i + k]* radius +clonedAtoms[0][4 * a + k]; //pos
@@ -200,15 +203,16 @@ if (!solid) {
 
 	delete[] tmp;
 	delete[] tmpi;
+
 	glBindVertexArray(0);
 	return e;
 } //SetupAtomsNoTess
 
 
-GLenum SetupAtoms(GLuint **AtomVAO, GLuint **AtomVertBuffer)
+GLenum SetupAtoms(GLuint **AtomVAO /*[3]*/, GLuint **AtomVertBuffer /*[2]*/, GLuint *BondIndices)
 {
 	if (!numAtoms)
-		return 0;
+		return glGetError();
 	//rgh FIXME: put this all in the same vao
 	
 	//http://prideout.net/blog/?p=48 //public domain code
@@ -221,10 +225,12 @@ GLenum SetupAtoms(GLuint **AtomVAO, GLuint **AtomVertBuffer)
 	}
 	eprintf("SetupAtoms: totalatoms=%d", totalatoms);
 
-	*AtomVAO = new GLuint[2]; //atoms, cloned atoms
+	*AtomVAO = new GLuint[3]; //atoms, cloned atoms, bonds
 	*AtomVertBuffer = new GLuint[2];
-	glGenVertexArrays(2, *AtomVAO);
+
+	glGenVertexArrays(3, *AtomVAO);
 	glGenBuffers(2, *AtomVertBuffer);
+	glGenBuffers(1, BondIndices);
 
 	glBindVertexArray((*AtomVAO)[0]);
 	glBindBuffer(GL_ARRAY_BUFFER, (*AtomVertBuffer)[0]);
@@ -233,26 +239,84 @@ GLenum SetupAtoms(GLuint **AtomVAO, GLuint **AtomVertBuffer)
 	glEnableVertexAttribArray(1);
 	glDisableVertexAttribArray(2);
 	glDisableVertexAttribArray(3);
+
+	e=glGetError();
+	if (e!=GL_NO_ERROR)
+		eprintf ("gl error %d, %s %d", e, __FILE__, __LINE__);
 	float *tmp = new float[4 * totalatoms];
 	float *current=tmp;
+	
+	const int atomlimit=30;
+
+	numBonds=new int[TIMESTEPS];
 	for (int p=0;p<TIMESTEPS;p++) {
 
-		for (int a = 0; a < numAtoms[p]; a++) {
-			for (int k = 0; k < 4; k++) {
-				*current++ = atoms[p][4 * a + k];
+			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]));
+						if (d*0.9f<r*r) {// bond
+							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++) {
+				m[k]=M[k]=tmp[k];
 			}
-		} //a
+			for (int a = 1; a < numAtoms[p]; a++) {
+				for (int k=0; k<3;k++) {
+					if (m[k]>tmp[4*a+k])
+						m[k]=tmp[4*a+k];
+					if (M[k]<tmp[4*a+k])
+						M[k]=tmp[4*a+k];
+				}
+			}
+			grid g(m, M, pow(numAtoms[p], 1.0/3), 0.9f);
+			for (int a = 1; a < numAtoms[p]; a++) 
+				g.add(tmp+4*a);
+			for (int a = 0; a < numAtoms[p]; a++) {
+				std::vector<float*> found=g.find(tmp+4*a);
+				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]));
+					bonds.push_back(((found[b]-tmp)/4)+(p==0?0:numAtoms[p-1]));
+				}
+			}
+		}
+		numBonds[p]=bonds.size();
 		if (p!=0)
 			numAtoms[p]+=numAtoms[p-1];
 	} //p
+
 	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,
-			GL_STATIC_DRAW);
-		if ((e = glGetError()) != GL_NO_ERROR)
-			eprintf( "opengl error %d, glBufferData, l %d\n", e, __LINE__);
+		GL_STATIC_DRAW);
+	if ((e = glGetError()) != GL_NO_ERROR)
+		eprintf( "opengl error %d, glBufferData, l %d\n", e, __LINE__);
 
-		glBindVertexArray(0);
+	glBindVertexArray(0);
 
 	if ((e = glGetError()) != GL_NO_ERROR)
 		eprintf( "opengl error %d, end of SetupAtoms, l %d\n", e, __LINE__);
@@ -273,22 +337,34 @@ GLenum SetupAtoms(GLuint **AtomVAO, GLuint **AtomVertBuffer)
 			atomtrajectoryrestarts[t].push_back(0);
 			for (int p=1;p<TIMESTEPS;p++) {
 				int a=atomtrajectories[t];
-				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);
+				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);
 			}
 			atomtrajectoryrestarts[t].push_back(TIMESTEPS);
 		}
 	}
 	delete[] tmp;
+	//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__);
 
 	//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(),
@@ -309,9 +385,11 @@ GLenum SetupAtoms(GLuint **AtomVAO, GLuint **AtomVertBuffer)
 
 GLenum SetupUnitCell(GLuint *UnitCellVAO, GLuint *UnitCellVertBuffer, GLuint *UnitCellIndexBuffer)
 {
-	if (!has_abc)
-		return 0;
 	GLenum e;
+	if ((e = glGetError()) != GL_NO_ERROR)
+		eprintf( "opengl error %d, begin of SetupUnitCell\n", e, __LINE__);
+	if (!has_abc)
+		return e;
 	glGenVertexArrays(1, UnitCellVAO);
 	glGenBuffers(1, UnitCellVertBuffer);
 	glGenBuffers(1, UnitCellIndexBuffer);
@@ -356,10 +434,11 @@ GLenum SetupUnitCell(GLuint *UnitCellVAO, GLuint *UnitCellVertBuffer, GLuint *Un
 	glBufferData(GL_ARRAY_BUFFER, sizeof(float) * 3*8 , tmp,
 			GL_STATIC_DRAW);
 	if ((e = glGetError()) != GL_NO_ERROR)
-		eprintf( "opengl error %d, glBufferData, l %d\n", e, __LINE__);
+		eprintf( "opengl error %d, glBufferData vertex, l %d\n", e, __LINE__);
 	glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 3 * sizeof(float), (const void *)(0));
 	glBufferData(GL_ELEMENT_ARRAY_BUFFER, sizeof(tmpi), tmpi, GL_STATIC_DRAW);
-
+	if ((e = glGetError()) != GL_NO_ERROR)
+		eprintf( "opengl error %d, glBufferData index, l %d\n", e, __LINE__);
 	return e;
 }
 

NOMADVRLib/atomsGL.h

@@ -4,7 +4,7 @@
 #include "MyGL.h"
 
 GLenum atomTexture(GLuint t);
-GLenum SetupAtoms(GLuint **AtomVAO, GLuint **AtomVertBuffer);
+GLenum SetupAtoms(GLuint **AtomVAO, GLuint **AtomVertBuffer, GLuint *BondIndices);
 GLenum SetupAtomsNoTess (GLuint **AtomVAO, GLuint **AtomVertBuffer, GLuint **AtomIndexBuffer);
 GLenum SetupUnitCell(GLuint *UnitCellVAO, GLuint *UnitCellVertBuffer, GLuint *UnitCellIndexBuffer);
 

OpenVR/TimestepData/hellovr_opengl.vcxproj

 <?xml version="1.0" encoding="utf-8"?>
-<Project DefaultTargets="Build" ToolsVersion="12.0" xmlns="http://schemas.microsoft.com/developer/msbuild/2003">
+<Project DefaultTargets="Build" ToolsVersion="14.0" xmlns="http://schemas.microsoft.com/developer/msbuild/2003">
   <ItemGroup Label="ProjectConfigurations">
     <ProjectConfiguration Include="Debug|Win32">
       <Configuration>Debug</Configuration>
@@ -29,7 +29,7 @@
     <ConfigurationType>Application</ConfigurationType>
     <UseDebugLibraries>true</UseDebugLibraries>
     <CharacterSet>Unicode</CharacterSet>
-    <PlatformToolset>v120</PlatformToolset>
+    <PlatformToolset>v140</PlatformToolset>
   </PropertyGroup>
   <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'" Label="Configuration">
     <ConfigurationType>Application</ConfigurationType>
@@ -42,7 +42,7 @@
     <UseDebugLibraries>false</UseDebugLibraries>
     <WholeProgramOptimization>true</WholeProgramOptimization>
     <CharacterSet>Unicode</CharacterSet>
-    <PlatformToolset>v120</PlatformToolset>
+    <PlatformToolset>v140</PlatformToolset>
   </PropertyGroup>
   <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'" Label="Configuration">
     <ConfigurationType>Application</ConfigurationType>
@@ -178,6 +178,7 @@
     <ClCompile Include="NOMADVRLib\atomsGL.cpp" />
     <ClCompile Include="NOMADVRLib\CompileGLShader.cpp" />
     <ClCompile Include="NOMADVRLib\ConfigFile.cpp" />
+    <ClCompile Include="NOMADVRLib\Grid.cpp" />
     <ClCompile Include="NOMADVRLib\IsoShaders.cpp" />
     <ClCompile Include="NOMADVRLib\IsosurfacesGL.cpp" />
     <ClCompile Include="NOMADVRLib\polyhedron.cpp" />
@@ -197,6 +198,7 @@
     <ClInclude Include="NOMADVRLib\CompileGLShader.h" />
     <ClInclude Include="NOMADVRLib\ConfigFile.h" />
     <ClInclude Include="NOMADVRLib\eprintf.h" />
+    <ClInclude Include="NOMADVRLib\Grid.h" />
     <ClInclude Include="NOMADVRLib\IsoShaders.h" />
     <ClInclude Include="NOMADVRLib\IsosurfacesGL.h" />
     <ClInclude Include="NOMADVRLib\MyGL.h" />

OpenVR/TimestepData/hellovr_opengl_main.cpp

@@ -216,6 +216,7 @@ private: // OpenGL bookkeeping
 	//for atoms
 	GLuint *m_glAtomVertBuffer; // [TIMESTEPS];
 	GLuint *m_unAtomVAO; //[TIMESTEPS];
+	GLuint BondIndices;
 
 	//for unit cells
 	GLuint m_glUnitCellVertBuffer; // primitive, possibly non-primitive in further. Deformed cube
@@ -327,6 +328,14 @@ void dprintf( const char *fmt, ... )
 //pure windows, no sdl
 void eprintf( const char *fmt, ... )
 {
+	static int numerrors=0;
+	numerrors++;
+	if (numerrors==25) {
+		MessageBoxA(0, "Max messages reached, no further reporting", "Warning", 0);
+	}
+	if (numerrors>15) {
+		return;
+	}
 	va_list args;
 	char buffer[ 2048 ];
 
@@ -619,11 +628,28 @@ bool CMainApplication::BInitGL()
 	if( !CreateAllShaders() )
 		return false;
 
+	GLenum e;
+
 	SetupTexturemaps();
+	e=glGetError();
+	if (e!=GL_NO_ERROR)
+		eprintf ("gl error %d, %s %d", e, __FILE__, __LINE__);
 	SetupScene();
+	e=glGetError();
+	if (e!=GL_NO_ERROR)
+		eprintf ("gl error %d, %s %d", e, __FILE__, __LINE__);
 	SetupCameras();
+	e=glGetError();
+	if (e!=GL_NO_ERROR)
+		eprintf ("gl error %d, %s %d", e, __FILE__, __LINE__);
 	SetupStereoRenderTargets();
+	e=glGetError();
+	if (e!=GL_NO_ERROR)
+		eprintf ("gl error %d, %s %d", e, __FILE__, __LINE__);
 	SetupDistortion();
+	e=glGetError();
+	if (e!=GL_NO_ERROR)
+		eprintf ("gl error %d, %s %d", e, __FILE__, __LINE__);
 	SetupDepthPeeling();
 
 	SetupRenderModels();
@@ -1311,7 +1337,7 @@ void CMainApplication::SetupUnitCell()
 void CMainApplication::SetupAtoms()
 {
 	GLenum e;
-	e=::SetupAtoms(&m_unAtomVAO, &m_glAtomVertBuffer);
+	e=::SetupAtoms(&m_unAtomVAO, &m_glAtomVertBuffer, &BondIndices);
 //	GLuint *vao, *buffer, *index;
 //	e=::SetupAtomsNoTess(&vao, &buffer, &index);
 	if (e!=GL_NO_ERROR)
@@ -2016,19 +2042,39 @@ if (numClonedAtoms!=0 && currentset==0) {
 		dprintf("Gl error after Render cloned Atom timestep =%d: %d, %s\n", currentset, e, gluErrorString(e));
 }
 
+//now bonds
+if (numBonds) {
+	glBindVertexArray(m_unAtomVAO[2]);
+	glUseProgram(m_unUnitCellProgramID);
+	glUniformMatrix4fv(m_nUnitCellMatrixLocation, 1, GL_FALSE, transform.get());
+	float color[4]={0.5,0.5,1,1};
+	glUniform4fv(m_nUnitCellColourLocation, 1, color);
+	if (currentset==0)
+		glDrawElements(GL_LINES, numBonds[0],  GL_UNSIGNED_INT, (void*)0);
+	else
+		glDrawElements(GL_LINES, numBonds[currentset]-numBonds[currentset-1], GL_UNSIGNED_INT, 
+			(void*)(sizeof(int)*numBonds[currentset-1]) );
+
+	if ((e = glGetError()) != GL_NO_ERROR)
+			dprintf("Gl error after Render Atom bonds timestep =%d: %d, %s\n", currentset, e, gluErrorString(e));
+}
+glBindVertexArray(m_unAtomVAO[0]);
+
 //now trajectories
 if (!showTrajectories)
 	return;
 
 glUseProgram(m_unUnitCellProgramID);
 glUniformMatrix4fv(m_nUnitCellMatrixLocation, 1, GL_FALSE, transform.get());
-float color[4]={1,0,0,1};
-glUniform4fv(m_nUnitCellColourLocation, 1, color);
+float color2[4]={1,0,0,1};
+glUniform4fv(m_nUnitCellColourLocation, 1, color2);
+
 if ((e = glGetError()) != GL_NO_ERROR)
 	dprintf("Gl error after glUniform4fv 2 RenderUnitCell: %d, %s\n", e, gluErrorString(e));
 glEnableVertexAttribArray(0);
 glDisableVertexAttribArray(1);
-
+if ((e = glGetError()) != GL_NO_ERROR)
+	dprintf("Gl error after Render Atom trajectories timestep =%d: %d, %s\n", currentset, e, gluErrorString(e));