From 2b38fe6af2d9aa85362a4626794e3b77914ddb30 Mon Sep 17 00:00:00 2001
From: Craig Fisher <66772456+cajfisher@users.noreply.github.com>
Date: Tue, 2 Nov 2021 18:19:49 +0900
Subject: [PATCH] Fix bugs in animation loop and arc file writer

Animation loop was returning to frame 1 instead of frame 0, and has now been fixed.
The arc writer also gave the wrong Cartesian coordinates for lattices rotated away from gdis's default, so matrix for calculating atom coordinates should be generated from arc's lattice parameters (pbcs).
---
 src/file_arc.c    | 370 ++++++++++++++++++++++++++--------------------
 src/gui_animate.c |   2 +-
 src/matrix.c      |  97 ++++++------
 src/matrix.h      |   3 +-
 4 files changed, 270 insertions(+), 202 deletions(-)

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 *);