diff --git a/src/file_arc.c b/src/file_arc.c
index 58641ef..c33c04f 100644
--- a/src/file_arc.c
+++ b/src/file_arc.c
@@ -45,7 +45,7 @@ extern struct elem_pak elements[];
 gint write_arc_frame(FILE *fp, struct model_pak *model)
 {
 gint m, flag;
-gdouble x[3];
+gdouble x[3], latmat[9];
 GSList *clist, *mlist;
 struct core_pak *core;
 struct mol_pak *mol;
@@ -53,14 +53,14 @@ time_t t1;
 
 /* get & print the date (streamlined slightly by sxm) */
 t1 = time(NULL);
-fprintf(fp,"!DATE %s",ctime(&t1));
+fprintf(fp, "!DATE %s", ctime(&t1));
 
 if (model->periodic)
   {
   if (model->periodic == 2)
     {
 /* saving a surface as a 3D model - make depth large enough to fit everything */
-    fprintf(fp,"PBC %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f (P1)\n",
+    fprintf(fp, "PBC %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f (P1)\n",
                 model->pbc[0], model->pbc[1], 2.0*model->rmax,
                 R2D*model->pbc[3],
                 R2D*model->pbc[4],
@@ -68,7 +68,7 @@ if (model->periodic)
     }
   else
     {
-    fprintf(fp,"PBC %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f (P1)\n",
+    fprintf(fp, "PBC %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f (P1)\n",
                 model->pbc[0], model->pbc[1], model->pbc[2],
                 R2D*model->pbc[3],
                 R2D*model->pbc[4],
@@ -77,14 +77,14 @@ if (model->periodic)
   }
 
 m=0;
-for (mlist=model->moles ; mlist ; mlist=g_slist_next(mlist))
+for (mlist = model->moles ; mlist ; mlist=g_slist_next(mlist))
   {
   mol = mlist->data;
 
 /* flag - atom written to file? (ie if UNDELETED) */
-  flag=0;
+  flag = 0;
 /* m - number of atoms in molecule */
-  for (clist=mol->cores ; clist ; clist=g_slist_next(clist))
+  for (clist = mol->cores ; clist ; clist=g_slist_next(clist))
     {
     core = clist->data;
     if (core->status & DELETED)
@@ -92,18 +92,20 @@ for (mlist=model->moles ; mlist ; mlist=g_slist_next(mlist))
 
 /* everything is cartesian after latmat mult */
     ARR3SET(x, core->x);
-    vecmat(model->latmat, x);
+/* use latmat of new lattice parameters to cope with rotated cells */
+    matrix_lattice_from_pbc(latmat, model->pbc, model->periodic);
+    vecmat(latmat, x);
 /* unique molecule numbering for BIOSYM files (>0) */
     if (core->atom_type)
       {
       fprintf(fp, "%-4s  %14.9f %14.9f %14.9f CORE %-6d %-7s %-2s ",
-                  core->atom_label,x[0],x[1],x[2],m+1, core->atom_type,
+                  core->atom_label, x[0], x[1], x[2], m+1, core->atom_type,
                   elements[core->atom_code].symbol);
       }
     else
       {
       fprintf(fp, "%-4s  %14.9f %14.9f %14.9f CORE %-6d ?       %-2s ",
-                  core->atom_label,x[0],x[1],x[2],m+1, 
+                  core->atom_label, x[0], x[1], x[2], m+1, 
                   elements[core->atom_code].symbol);
       }
 
@@ -118,11 +120,11 @@ for (mlist=model->moles ; mlist ; mlist=g_slist_next(mlist))
 /* if at least one atom in the molecule was written - inc mol numbering */
   if (flag)
     {
-    fprintf(fp,"end\n");
+    fprintf(fp, "end\n");
     m++;
     }
   }
-fprintf(fp,"end\n");
+fprintf(fp, "end\n");
 
 return(0);
 }
@@ -133,11 +135,11 @@ return(0);
 gint write_arc_header(FILE *fp, struct model_pak *model)
 {
 /* print header */
-fprintf(fp,"!BIOSYM archive 3\n");
+fprintf(fp, "!BIOSYM archive 3\n");
 if (model->periodic)
-  fprintf(fp,"PBC=ON\n");
+  fprintf(fp, "PBC=ON\n");
 else
-  fprintf(fp,"PBC=OFF\n");
+  fprintf(fp, "PBC=OFF\n");
 gdis_blurb(fp);
 
 return(0);
@@ -213,8 +215,8 @@ struct shel_pak *shel;
 gchar *line,*tamp;
 gdouble energy;
 #define DEBUG_ARC 0
-gchar c_type[5],c_ptype[8],c_symb[3];
-gdouble x,y,z,charge;
+gchar c_type[5], c_ptype[8], c_symb[3];
+gdouble x, y, z, charge;
 
 g_assert(fp != NULL);
 g_assert(data != NULL);
@@ -229,154 +231,202 @@ clist = data->cores;
 slist = data->shels;
 end_count=0;
 
-	/*The energy is stored on the frame 1st line starting from the 65th caracter*/
+/* energy is stored on the 1st line of frame from 65th charactear */
 line = file_read_line(fp);
 sscanf(&(line[65])," %lf",&energy);
 //this now trigger a heap-buffer-overflow, replacing with safer g_strdup_printf version
-//sprintf(line,"%lf eV",energy);/*is it really always eV?*/
-//g_sprintf(line,"%lf eV",energy);
-//property_add_ranked(3,"Energy",line,data);
-tamp=g_strdup_printf("%lf eV",energy);
-property_add_ranked(3,"Energy",tamp,data);
+//sprintf(line,"%lf eV", energy); /* is it really always eV? */
+//g_sprintf(line,"%lf eV", energy);
+//property_add_ranked(3, "Energy", line, data);
+tamp=g_strdup_printf("%lf eV", energy);
+property_add_ranked(3, "Energy", tamp, data);
 g_free(tamp);
-	/*Note that title of calculation is given in the first 64 character BUT it can be empty*/
-/*get to the begining*/
-while(g_ascii_strncasecmp("!date", line, 5)!=0) {
-	g_free(line);
-	line = file_read_line(fp);
-}
-g_free(line);/*should we do something with the !DATE line?*/
+/* Nb. title of calculation is given in the first 64 characters BUT it can be empty*/
+/* move to the begining */
+while(g_ascii_strncasecmp("!date", line, 5) != 0)
+  {
+  g_free(line);
+  line = file_read_line(fp);
+  }
+g_free(line); /* should we do something with the !DATE line? */
 /* loop while there's data */
-for (;;) {
-/*we do not use the tokenized method to avoid some break in threaded environment*/
+for (;;)
+  {
+/* we do not use the tokenized method to avoid some break in threaded environment */
   line = file_read_line(fp);
-  if(!line) break;
+  if(!line)
+    break;
 
-/* increment/reset end counter according to input line <- strstr at the begining of the line*/
-  if (g_ascii_strncasecmp("end", line, 3) == 0) end_count++;
-  else end_count = 0;
+/* increment/reset end counter according to input line <- strstr at the begining of the line */
+  if (g_ascii_strncasecmp("end", line, 3) == 0)
+    end_count++;
+  else
+    end_count = 0;
 /* skip single end, terminate on double end */
-  if (end_count == 1) {
-	g_free(line);
-	continue;
-  }
-  if (end_count == 2) break;/*normal exit*/
+  if (end_count == 1)
+    {
+    g_free(line);
+    continue;
+    }
+  if (end_count == 2)
+    break; /* normal exit */
 
 /* cell dimension search */
-  if (g_ascii_strncasecmp("pbc", line, 3) == 0){
+  if (g_ascii_strncasecmp("pbc", line, 3) == 0)
+    {
     data->pbc[5]=-1.;
-sscanf(line,"PBC %lf %lf %lf %lf %lf %lf %*s",&(data->pbc[0]),&(data->pbc[1]),&(data->pbc[2]),&(data->pbc[3]),&(data->pbc[4]),&(data->pbc[5]));
+    sscanf(line, "PBC %lf %lf %lf %lf %lf %lf %*s", &(data->pbc[0]), &(data->pbc[1]), 
+		    &(data->pbc[2]), &(data->pbc[3]), &(data->pbc[4]), &(data->pbc[5]));
 #if DEBUG_ARC
-fprintf(stdout,"#DBG: LINE:  %lf %lf %lf %lf %lf %lf\n",data->pbc[0],data->pbc[1],data->pbc[2],data->pbc[3],data->pbc[4],data->pbc[5]);
+fprintf(stdout,"#DBG: LINE:  %lf %lf %lf %lf %lf %lf\n", data->pbc[0], data->pbc[1], 
+		             data->pbc[2], data->pbc[3], data->pbc[4], data->pbc[5]);
 #endif
-    if((data->pbc[5])>=0){/*6 tokens -> PCB=ON with or without implicit 2D*/
-		data->pbc[3] *= D2R;
-		data->pbc[4] *= D2R;
-		data->pbc[5] *= D2R;
-		if((data->pbc[2])==0){/*implicit 2D*/
-			data->pbc[2] = 1.0;
-			data->pbc[3] = 0.5*G_PI;
-			data->pbc[4] = 0.5*G_PI;
-			data->periodic = 2;
-		}
-    }else{/*less than 6 tokens -> PCB=2D */
-		data->pbc[5] = D2R*(data->pbc[2]);
-		data->pbc[2] = 1.0;
-		data->pbc[3] = 0.5*G_PI;
-		data->pbc[4] = 0.5*G_PI;
-    }
+    if((data->pbc[5])>=0)
+      { /* 6 tokens -> PCB=ON with or without implicit 2D */
+	data->pbc[3] *= D2R;
+	data->pbc[4] *= D2R;
+	data->pbc[5] *= D2R;
+	if((data->pbc[2])==0)
+          { /* implicit 2D */
+	  data->pbc[2] = 1.0;
+	  data->pbc[3] = 0.5*G_PI;
+	  data->pbc[4] = 0.5*G_PI;
+	  data->periodic = 2;
+	  }
+        }
+      else
+        { /* less than 6 tokens -> PCB=2D */
+	data->pbc[5] = D2R*(data->pbc[2]);
+	data->pbc[2] = 1.0;
+	data->pbc[3] = 0.5*G_PI;
+	data->pbc[4] = 0.5*G_PI;
+        }
+      matrix_lattice_init(data); /* initialize matrix */
 #if DEBUG_ARC
-fprintf(stdout,"#DBG: a=%lf b=%lf c=%lf alpha=%lf beta=%lf gamma=%lf\n",data->pbc[0],data->pbc[1],data->pbc[2],data->pbc[3],data->pbc[4],data->pbc[5]);
+fprintf(stdout, "#DBG: a=%lf b=%lf c=%lf alpha=%lf beta=%lf gamma=%lf\n", 
+	data->pbc[0], data->pbc[1], data->pbc[2], data->pbc[3], data->pbc[4], data->pbc[5]);
+      P3MAT("lattice matrix", model->latmat);
 #endif
-    matrix_lattice_init(data);/*initialize matrix now*/
-    data->fractional = TRUE;/*set fractional now*/
-    g_free(line);
-    continue;
-  }
+      data->fractional = TRUE; /* set fractional */
+      g_free(line);
+      continue;
+    }
 
-/* each coord line is defined as:
+/* each coord line consists of
  * [name], x, y, z, type, [seq. nb], pot. type, elem symbol, charge [, atom nb] */
-  c_symb[0]='-';c_symb[1]='\0';
-  sscanf(line,"%*s %lf %lf %lf %s %*s %s %s %lf %*s",&x,&y,&z,&(c_type[0]),&(c_ptype[0]),&(c_symb[0]),&charge);
+    c_symb[0]='-';c_symb[1]='\0';
+    sscanf(line,"%*s %lf %lf %lf %s %*s %s %s %lf %*s", 
+		  &x, &y, &z, &(c_type[0]), &(c_ptype[0]), &(c_symb[0]), &charge);
 #if DEBUG_ARC
-fprintf(stdout,"#DBG: LINE: %lf %lf %lf %s %s %s %f\n",x,y,z,c_type,c_ptype,c_symb,charge);
+fprintf(stdout,"#DBG: LINE: %lf %lf %lf %s %s %s %f\n", x, y, z, c_type, c_ptype, c_symb, charge);
 #endif
-  if(c_symb[0]=='-'){
-	gui_text_show(ERROR, "unexpected end of file reading arc file\n");
-	g_free(line);
-	return;
-  }
+    if(c_symb[0]=='-')
+      {
+      gui_text_show(ERROR, "Unexpected end of file reading arc file\n");
+      g_free(line);
+      return;
+      }
 /* process core/shell candidates */
-  if (elem_symbol_test(c_symb)){
-	core_flag=TRUE;
-	region = 0;
+    if (elem_symbol_test(c_symb))
+      {
+      core_flag=TRUE;
+      region = 0;
 /* MARVIN core/shell/region parse */
-	if (is_marvin_label(c_type)){
-		if (!marvin_core(c_type)) core_flag = FALSE;
-		region = marvin_region(c_type);
-		data->region_empty[region] = FALSE;
-	}
+      if (is_marvin_label(c_type))
+        {
+        if (!marvin_core(c_type)) core_flag = FALSE;
+          region = marvin_region(c_type);
+        data->region_empty[region] = FALSE;
+        }
 /* overwrite existing core (if any) */
-	if (core_flag){
-		if (clist){
-			core = clist->data;
-			clist = g_slist_next(clist);
-		}else{
+      if (core_flag)
+        {
+        if (clist){
+	  core = clist->data;
+	  clist = g_slist_next(clist);
+        }
+      else
+        {
 /* otherwise, add a new core */
 /* NB: must append (not prepend) since we may overwrite & the order must be the same */
 /* TODO - we could prepend cores (it's faster) and then reverse the list at the end */
 /* but it's more complicated as we should ONLY reverse newly added cores */
-			core=core_new(c_symb,"---",data);/*to make sure we don't need to free/dup atom_label all the time*/
-			core->atom_type=g_malloc(9*sizeof(gchar));/*same here for atom_type*/
-			data->cores = g_slist_append(data->cores, core);
-		}
+        core=core_new(c_symb, "---", data); /* to make sure we don't need to free/dup atom_label all the time */
+        core->atom_type=g_malloc(9*sizeof(gchar)); /* same here for atom_type */
+        data->cores = g_slist_append(data->cores, core);
+        }
 	
 /* set the proper label */
-	sprintf(core->atom_label,"%s",c_symb);/*FIX 7637ff*/
+	sprintf(core->atom_label, "%s", c_symb); /* FIX 7637ff */
 
 /* set values */
-	core->x[0] = x;core->x[1] = y;core->x[2] = z;core->charge = charge;
-if(data->periodic) {
-	/* convert input cartesian coords to fractional now */
-	vecmat(data->ilatmat, core->x);
-	if(core->x[0]<0) core->x[0]+=1.0;
-	if(core->x[1]<0) core->x[1]+=1.0;
-	if(core->x[2]<0) core->x[2]+=1.0;
-	if(core->x[0]>1.0) core->x[0]-=1.0;
-	if(core->x[1]>1.0) core->x[1]-=1.0;
-	if(core->x[2]>1.0) core->x[2]-=1.0;
-}
-	sprintf(core->atom_type,"%s",c_ptype);
-	core->region = region;core->lookup_charge = FALSE;
-	}/*FIX 7637ff (ugly one: in case core_flag==FALSE core=NULL and all core->xxx FAIL)*/	
-  }else{
+	core->x[0] = x;
+	core->x[1] = y;
+	core->x[2] = z;
+	core->charge = charge;
+    if(data->periodic)
+      {
+      /* convert input cartesian coords to fractional now */
+      vecmat(data->ilatmat, core->x);
+      if(core->x[0] < 0)
+        core->x[0] += 1.0;
+      if(core->x[1] < 0)
+        core->x[1] += 1.0;
+      if(core->x[2] < 0)
+        core->x[2] += 1.0;
+      if(core->x[0] > 1.0)
+        core->x[0] -= 1.0;
+      if(core->x[1] > 1.0)
+        core->x[1] -= 1.0;
+      if(core->x[2] > 1.0)
+        core->x[2] -= 1.0;
+      }
+    sprintf(core->atom_type,"%s",c_ptype);
+    core->region = region;core->lookup_charge = FALSE;
+    } /* FIX 7637ff (ugly one: in case core_flag==FALSE core=NULL and all core->xxx FAIL) */	
+  }
+else
+  {
 /* overwrite existing shell (if any) */
-	if (slist){
-        	shel = slist->data;
-        	slist = g_slist_next(slist);
-        }else{
+  if (slist)
+    {
+    shel = slist->data;
+    slist = g_slist_next(slist);
+    }
+  else
+    {
 /* otherwise, add a new shell */
-		shel=shell_new(c_symb,"---",data);/*to make sure we don't need to free/dup shell_label all the time*/
-		data->shels = g_slist_append(data->shels, shel);
-	}
+    shel = shell_new(c_symb, "---", data); /* to make sure we don't need to free/dup shell_label all the time */
+    data->shels = g_slist_append(data->shels, shel);
+    }
 /* set the proper label */
-  sprintf(shel->shell_label,"%s",c_symb);
+    sprintf(shel->shell_label, "%s", c_symb);
 /* set values */
-  shel->x[0] = x;shel->x[1] = y;shel->x[2] = z;shel->charge = charge;
-if(data->periodic) {
-	/* convert input cartesian coords to fractional now */
-	vecmat(data->ilatmat, shel->x);
-	if(shel->x[0]<0) shel->x[0]+=1.0;
-	if(shel->x[1]<0) shel->x[1]+=1.0;
-	if(shel->x[2]<0) shel->x[2]+=1.0;
-	if(shel->x[0]>1.0) shel->x[0]-=1.0;
-	if(shel->x[1]>1.0) shel->x[1]-=1.0;
-	if(shel->x[2]>1.0) shel->x[2]-=1.0;
-}
-  shel->region = region;shel->lookup_charge = FALSE;
+    shel->x[0] = x;
+    shel->x[1] = y;
+    shel->x[2] = z;
+    shel->charge = charge;
+    if(data->periodic)
+      {
+      /* convert input cartesian coords to fractional now */
+      vecmat(data->ilatmat, shel->x);
+      if(shel->x[0] < 0)
+        shel->x[0] += 1.0;
+      if(shel->x[1] < 0)
+        shel->x[1] += 1.0;
+      if(shel->x[2] < 0)
+        shel->x[2] += 1.0;
+      if(shel->x[0] > 1.0)
+        shel->x[0] -= 1.0;
+      if(shel->x[1] > 1.0)
+        shel->x[1] -= 1.0;
+      if(shel->x[2] > 1.0)
+        shel->x[2] -= 1.0;
+      }
+    shel->region = region;shel->lookup_charge = FALSE;
+    }
+    g_free(line);
   }
-  g_free(line);
-}
 g_free(line);
 }
 
@@ -413,47 +463,53 @@ if (!fp)
 /* loop while there's data */
 for (;;)
   {
-  fp_pos=ftell(fp);/*get the past line*/
+  fp_pos = ftell(fp); /* get position of line before reading next */
   line = file_read_line(fp);
-  if(!line) break;
+  if(!line)
+    break;
 
 /* pbc on/off search */
-if(g_ascii_strncasecmp("pbc=on", line, 6) == 0) {
-	data->periodic = 3;
-	old_fp_pos=fp_pos;
-	g_free(line);
-	continue;
-}
-if(g_ascii_strncasecmp("pbc=off", line, 7) == 0) {
-	data->periodic = 0;
-	old_fp_pos=fp_pos;
-	g_free(line);
-	continue;
-}
-if(g_ascii_strncasecmp("pbc=2d", line, 6) == 0) {
-	data->periodic = 2;
-	old_fp_pos=fp_pos;
-	g_free(line);
-	continue;
-}
-/* coords search */
-if(g_ascii_strncasecmp("!date", line, 5) == 0) {
+  if(g_ascii_strncasecmp("pbc=on", line, 6) == 0)
+    {
+    data->periodic = 3;
+    old_fp_pos = fp_pos;
+    g_free(line);
+    continue;
+    }
+  if(g_ascii_strncasecmp("pbc=off", line, 7) == 0)
+    {
+    data->periodic = 0;
+    old_fp_pos = fp_pos;
+    g_free(line);
+    continue;
+    }
+  if(g_ascii_strncasecmp("pbc=2d", line, 6) == 0)
+    {
+    data->periodic = 2;
+    old_fp_pos = fp_pos;
+    g_free(line);
+    continue;
+    }
+  /* coords search */
+  if(g_ascii_strncasecmp("!date", line, 5) == 0)
+    {
 /* NEW */
-    fp_pos=ftell(fp);/*tag here*/
-    fseek(fp,old_fp_pos,SEEK_SET);/* rewind a line before !date (to get energy) */
+    fp_pos = ftell(fp); /* tag here */
+    fseek(fp, old_fp_pos, SEEK_SET); /* rewind a line before !date (to get energy) */
     add_frame_offset(fp, data);
     
 /* only load the required frame */
     if (frame == data->cur_frame)
       read_arc_block(fp, data);
-    fseek(fp,fp_pos,SEEK_SET);/*go just after !date*/
+    fseek(fp, fp_pos, SEEK_SET); /* go just after !date */
 /* increment counter */
     frame++;
     }
   g_free(line);
-  old_fp_pos=fp_pos;
+  old_fp_pos = fp_pos;
   }
-g_free(line);
+//g_free(line);
+
 /* got everything */
 data->num_frames = frame;
 
diff --git a/src/gui_animate.c b/src/gui_animate.c
index 970ca89..3b9f4a0 100644
--- a/src/gui_animate.c
+++ b/src/gui_animate.c
@@ -256,7 +256,7 @@ g_assert(model != NULL);
 model->cur_frame += model->anim_step;
 if (model->cur_frame >= model->num_frames)
   if (model->anim_loop)
-    model->cur_frame = 1;
+    model->cur_frame = 0;
 
 /* continue until we run out of frames (or a stop is flagged) */
 if (model->cur_frame < model->num_frames && model->animating)
diff --git a/src/matrix.c b/src/matrix.c
index 61bc73c..c5d3c16 100644
--- a/src/matrix.c
+++ b/src/matrix.c
@@ -42,9 +42,9 @@ The GNU GPL can also be found at http://www.gnu.org
 /* main structure */
 extern struct sysenv_pak sysenv;
 
-/**********************************************************/
-/* adjust a rotation matrix wrt curent camera position */
-/**********************************************************/
+/********************************************************/
+/* adjust a rotation matrix wrt current camera position */
+/********************************************************/
 void matrix_camera_transform(gdouble *rot, struct model_pak *model)
 {
 gdouble q[9], qi[9], mat[9];
@@ -752,6 +752,53 @@ switch(type)
   }
 }
 
+/**********************************************/
+/* generate standard lattice matrix from pbcs */
+/**********************************************/
+void matrix_lattice_from_pbc(gdouble *mat, gdouble *pbc, gint periodic)
+{
+gdouble b1, b2, b3, c1, c2, c3;
+
+g_assert(mat != NULL);
+g_assert(pbc != NULL);
+
+/* FIXME - correction needed for 1D case as well? */
+if (periodic == 2)
+  {
+/* lattice code requires non-existent c parameter to be 1 */
+  pbc[2] = 1.0;
+  }
+
+/* compute the translation vector for b */
+b1 = cos(pbc[5]);
+b2 = sin(pbc[5]);
+b3 = 0.0;         /* constrain a,b to remain on the x-y plane */
+
+if (b2 == 0.0)
+  {
+  printf("matrix_lattice_from_pbc(): bad cell parameters.\n");
+  b2 = 1.0;
+  }
+
+/* compute the translation vector for c */
+c1 = cos(pbc[4]);
+c2 = (2.0*cos(pbc[3]) + b1*b1 + b2*b2 - 2.0*b1*c1 - 1.0)/(2.0*b2);
+c3 = sqrt(1.0 - c1*c1 - c2*c2);
+
+/* assign in rows to make it easier to scale */
+/* x & a are assumed co-linear */
+VEC3SET(&mat[0], pbc[0], 0.0, 0.0);
+VEC3SET(&mat[3], b1, b2, b3);
+VEC3SET(&mat[6], c1, c2, c3);
+
+/* scale b & c vectors up */
+VEC3MUL(&mat[3], pbc[1]);
+VEC3MUL(&mat[6], pbc[2]);
+
+/* get vectors in cols */
+matrix_transpose(mat);
+}
+
 /**************************************************************/
 /* get reciprocal lattice vectors from direct lattice vectors */
 /**************************************************************/
@@ -803,7 +850,7 @@ data->rlatmat[5] = c[1];
 data->rlatmat[8] = c[2];
 
 #if DEBUG_RLAT
-P3MAT("reciprocal lattice matrix:",data->rlatmat);
+P3MAT("reciprocal lattice matrix:", data->rlatmat);
 #endif
 }
 
@@ -814,10 +861,9 @@ P3MAT("reciprocal lattice matrix:",data->rlatmat);
 void matrix_lattice_init(struct model_pak *data)
 {
 gdouble n[3], v1[3], v2[3], v3[3];
-gdouble b1, b2, b3, c1, c2, c3;
 
 /* use a supplied latmat, rather than generating from pbc's */
-/* NB: should be in gdis style colum vector format (gulp/marvin are in rows) */
+/* NB: should be in gdis style column vector format (gulp/marvin are in rows) */
 if (data->construct_pbc)
   {
 #if DEBUG_XLAT
@@ -855,47 +901,14 @@ printf("constructing pbc...\n");
 else
   {
 #if DEBUG_XLAT
-printf("constructing latmat...\n");
+  printf("constructing latmat...\n");
 #endif
 /* construct lattice matrix from the unit cell lengths & angles */
 /* this basically works by using the cosine rule in conjunction */
 /* with a few geometric constraints (eg normalized vectors) */
 
 /* FIXME - correction needed for 1D case as well? */
-if (data->periodic == 2)
-  {
-/* lattice code requires non existent c parameter to be 1 */
-  data->pbc[2] = 1.0;
-  }
-
-/* compute the translation vector for b */
-  b1 = cos(data->pbc[5]);
-  b2 = sin(data->pbc[5]);
-  b3 = 0.0;              /* constrain a,b to remain on the x,y plane */
-
-  if (b2 == 0.0)
-    {
-    printf("matrix_lattice_init(): bad cell parameters.\n");
-    b2 = 1.0;
-    }
-
-/* compute the translation vector for c */
-  c1 = cos(data->pbc[4]);
-  c2 = (2.0*cos(data->pbc[3]) + b1*b1 + b2*b2 - 2.0*b1*c1 - 1.0)/(2.0*b2);
-  c3 = sqrt(1.0 - c1*c1 - c2*c2);
-
-/* assign in rows to make it easier to scale */
-/* x & a are assumed co-linear */
-  VEC3SET(&data->latmat[0], data->pbc[0], 0.0, 0.0);
-  VEC3SET(&data->latmat[3], b1, b2, b3);
-  VEC3SET(&data->latmat[6], c1, c2, c3);
-
-/* scale b & c vectors up */
-  VEC3MUL(&data->latmat[3], data->pbc[1]);
-  VEC3MUL(&data->latmat[6], data->pbc[2]);
-
-/* get vectors in cols */
-  matrix_transpose(data->latmat);
+  matrix_lattice_from_pbc(data->latmat, data->pbc, data->periodic);
   }
 
 /* update dependents */
@@ -922,7 +935,6 @@ VEC3MUL(data->cell_angles, R2D);
 printf("cell dimensions: [%6.2f %6.2f %6.2f] (%6.2f %6.2f %6.2f)\n",
        data->pbc[0], data->pbc[1], data->pbc[2],
        data->cell_angles[0], data->cell_angles[1], data->cell_angles[2]);
-P3MAT("lattice matrix:",data->latmat);
 #endif
 return;
 }
@@ -950,7 +962,6 @@ g_assert(model != NULL);
 
 #if DEBUG_NEW_LATTICE
 P3MAT("transformation matrix: :", latmat);
-P3MAT("old latmat:", model->latmat);
 #endif
 
 /* get the transformation matrix */
diff --git a/src/matrix.h b/src/matrix.h
index c9dd9b2..41befdb 100644
--- a/src/matrix.h
+++ b/src/matrix.h
@@ -52,10 +52,11 @@ void matrix_transpose_44(gdouble *);
 
 void calc_norm(gdouble *, gdouble *, gdouble *, gdouble *);
 
+void matrix_lattice_from_pbc(gdouble *, gdouble *, gint);
 void matrix_lattice_init(struct model_pak *);
 void matrix_lattice_new(gdouble *, struct model_pak *);
 void matrix_identity(gdouble *);
-void matrix_rotation(gdouble *, gdouble , gint);
+void matrix_rotation(gdouble *, gdouble, gint);
 
 void matrix_z_alignment(gdouble *, gdouble *);
 void matrix_v_alignment(gdouble *, gdouble *, gdouble *);