diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md
index 727a9cb..4fe1c87 100644
--- a/.github/CHANGELOG.md
+++ b/.github/CHANGELOG.md
@@ -3,6 +3,7 @@
 ### New features
 
 * Add partial extended xyz support (https://github.com/briling/v/pull/18)
+* Read and print lattice parameters to/from extxyz header (#39)
 * Add `j` hotkey to jump to a frame inside and CLI option `frame:%d` to start with a specific frame (https://github.com/briling/v/pull/5)
 * Add `bmax` CLI argument for max. bond length (920202b)
 * Add `com` and `exitcom` CLI argument for `gui:0` and on-exit command sequences, respectively
@@ -17,6 +18,7 @@
 * Read molecules from the standard input (#28)
 * Show infrared intensities and mode masses (#35)
 * Improve the text-in-corner look (de242c7, #36)
+* Add `u` to the headless mode (#37)
 * Improve data structures and code readability
 
 ### Fixes
@@ -28,7 +30,6 @@
 * Fix text blinking when playing animation (80cae88, #36)
 
 ### Coming in the next version:
-* extended xyz (#16, #17)
 * high-symmetry determination bugs (#21)
 * how to build on mac (#34)
 
diff --git a/README.md b/README.md
index 525a1fd..b546187 100644
--- a/README.md
+++ b/README.md
@@ -45,13 +45,14 @@ To build:
 ```
 ./v file [file2 ... fileN] [options]
 ```
-A filename `-` stands for the standard input.
+A filename `-` stands for the standard input (xyz files only).
 
 Show the reference:
 ```
 ./v
 ```
 
+### Options
 <details open><summary><strong>Command-line options</strong></summary>
 
 |                          |                                                               |
@@ -81,6 +82,7 @@ Show the reference:
 
 </details>
 
+### Keyboard
 <details open><summary><strong>Keyboard reference</strong></summary>
 
 |                                |                                                           |
@@ -120,15 +122,17 @@ Show the reference:
 
 </details>
 
-<details open><summary><strong>Mouse (in development)</strong></summary>
+### Mouse
+
 One can also use the mouse to rotate the molecule and zoom in/out.
-</details>
 
-<details open><summary><strong>Headless mode (in development)</strong></summary>
+### Additional commands
 
-If run in the headless mode with `gui:0`, the symbols from the standard input are processed 
+#### Headless mode
+
+If run in the headless mode with `gui:0`, the symbols from the standard input are processed
 as if the corresponding keys were pressed in the normal mode.
-Right now, only `p`, `x`, `z`, and `.` are available.
+Right now, `p`, `x`, `z`, `u`, and `.` are available.
 Command-line option `com:%s` overrides the standard input.
 These examples are equivalent:
 ```
@@ -142,8 +146,60 @@ D*h
 D*h
 ```
 
+#### Normal mode
+In the normal mode, the symbols from the CLI option `exitcom:` are executed immediately before closing.
+For example,
+```
+./v mol/mol0001.xyz exitcom:z
+```
+automatically prints the last xyz coordinates when the user closes the window.
+
+
+
+
+### Boundary conditions
+Two types of boundary conditions are recognized:
+* PBC (3D)
+* spherical confinement
+
+#### PBC
+PBC can be read from the xyz [file header](https://ase-lib.org/ase/io/formatoptions.html#extxyz)
+by specifying `Lattice="ax ay az bx by bz cx cy cz"`:
+```
+./v mol/MOL_3525.ext.xyz
+```
+Currently, only the PBC in all three dimensions are supported.
+
+The same can be passed via the command-line, in which case it overrides the one from the file:
+```
+./v mol/MOL_3525.xyz cell:8.929542,0,0,4.197206,8.892922,0,0.480945,2.324788,10.016044
+```
+For orthogonal/cubic cell:
+```
+./v mol/1372_D02.340_1.out bonds:0 cell:20.23,20.23,20.23
+./v mol/1372_D02.340_1.out bonds:0 cell:20.23
+```
+In Bohr instead of Å:
+```
+./v mol/1372_D02.340_1.out bonds:0 cell:b10.7
+```
+Finally, to disable the cell from the file:
+```
+./v mol/MOL_3525.ext.xyz cell:0
+```
+
+#### Spherical confinement
+Spherical confinement can be specified from the command-line by the following:
+```
+./v mol/mol0001.xyz shell:2     #  sphere with r = 2 Å is put around the molecule
+./v mol/mol0001.xyz shell:b4    #  sphere with r = 2 Bohr
+./v mol/mol0001.xyz shell:2,3   #  spheres with r = 2 and 3 Å (e.g., soft and hard boundaries)
+./v mol/mol0001.xyz shell:b4,5  #  spheres with r = 4 and 5 Bohr
+```
+
 </details>
 
+
 ## Examples [↑](#contents)
 * `mol/C3H6~mCPBA_01x11.qm.out` — geometries + vibrations
 ```
@@ -166,7 +222,7 @@ D*h
 ![Adamantane mode animation](figures/C10H16.qm.out.gif)
 * `mol/1372_D02.340_1.out` — PBC simulation
 ```
-./v mol/1372_D02.340_1.out bonds:0 cell:b10.7,10.7,1.07 font:-*-*-medium-*-*--15-*-*-*-*-*-*-1
+./v mol/1372_D02.340_1.out bonds:0 cell:b10.7 font:-*-*-medium-*-*--15-*-*-*-*-*-*-1
 ```
 ![Atoms in cell with PBC](figures/1372_D02.340_1.out_1087.gif)
 * `mol/mol0001.xyz`, `mol/mol0002.xyz` — `.xyz` files with atomic numbers and atomic symbols
@@ -177,6 +233,9 @@ D*h
 
 * `mol/MOL_3525.xyz` — organic crystal with non-orthogonal cell
 ```
+./v mol/MOL_3525.ext.xyz
+```
+```
 ./v mol/MOL_3525.xyz cell:8.929542,0.0,0.0,4.197206,8.892922,0.0,0.480945,2.324788,10.016044 font:-*-*-medium-*-*--15-*-*-*-*-*-*-1
 ```
 ![Organic crystal cell](figures/MOL_3525.xyz_1.gif)
diff --git a/obj/v/ac3_print.d b/obj/v/ac3_print.d
index 0ca5016..78de98c 100644
--- a/obj/v/ac3_print.d
+++ b/obj/v/ac3_print.d
@@ -1,2 +1,3 @@
 obj/v/ac3_print.o obj-pic/v/ac3_print.o: src/v/ac3_print.c src/v/v.h \
- src/mol/mol.h src/mol/common.h src/v/pars.h src/math/vec3.h
+ src/mol/mol.h src/mol/common.h src/v/pars.h src/math/matrix.h \
+ src/math/vecn.h src/math/vec3.h
diff --git a/obj/v/ac3_read.d b/obj/v/ac3_read.d
index 08c3d55..3b18538 100644
--- a/obj/v/ac3_read.d
+++ b/obj/v/ac3_read.d
@@ -1,3 +1,3 @@
 obj/v/ac3_read.o obj-pic/v/ac3_read.o: src/v/ac3_read.c src/v/v.h \
  src/mol/mol.h src/mol/common.h src/v/pars.h src/math/vecn.h \
- src/math/vec3.h
+ src/math/vec3.h src/math/matrix.h src/math/vecn.h
diff --git a/obj/v/ac3_read_xyz.d b/obj/v/ac3_read_xyz.d
index 9faf7c0..ccd0291 100644
--- a/obj/v/ac3_read_xyz.d
+++ b/obj/v/ac3_read_xyz.d
@@ -1,2 +1,2 @@
 obj/v/ac3_read_xyz.o obj-pic/v/ac3_read_xyz.o: src/v/ac3_read_xyz.c \
- src/v/v.h src/mol/mol.h src/mol/common.h src/v/pars.h
+ src/v/v.h src/mol/mol.h src/mol/common.h src/v/pars.h src/math/vecn.h
diff --git a/obj/v/cli.d b/obj/v/cli.d
index b83ea3a..25c9f99 100644
--- a/obj/v/cli.d
+++ b/obj/v/cli.d
@@ -1,3 +1,3 @@
 obj/v/cli.o obj-pic/v/cli.o: src/v/cli.c src/v/v.h src/mol/mol.h \
  src/mol/common.h src/v/pars.h src/math/vecn.h src/math/matrix.h \
- src/math/vecn.h src/math/vec3.h
+ src/math/vecn.h
diff --git a/src/v/ac3_draw.c b/src/v/ac3_draw.c
index d06b405..664b166 100644
--- a/src/v/ac3_draw.c
+++ b/src/v/ac3_draw.c
@@ -28,20 +28,20 @@ static int cmpz(const void * p1, const void * p2){
     return  0;
 }
 
-void ac3_draw(atcoord * ac, rendpars rend){
+void ac3_draw(atcoord * ac, rendpars * rend){
   int n = ac->n;
   kzstr * kz = malloc(sizeof(kzstr)*n);
-  int   * ks = (rend.bonds>0) ? malloc(sizeof(int)*n) : NULL;
+  int   * ks = (rend->bonds>0) ? malloc(sizeof(int)*n) : NULL;
 
   double resol = world.size * RESOL_SCALE;
-  double r1  = rend.r * resol * rend.scale;
+  double r1  = rend->r * resol * rend->scale;
 
   for(int k=0; k<n; k++){
     kz[k].k = k;
     kz[k].z = ac->r[k*3+2];
   }
   qsort(kz, n, sizeof(kzstr), cmpz);
-  if(rend.bonds>0){
+  if(rend->bonds>0){
     for(int i=0; i<n; i++){
       ks[ kz[i].k ] = i;
     }
@@ -61,19 +61,19 @@ void ac3_draw(atcoord * ac, rendpars rend){
       XDrawArc(world.dis, world.canv, q>0?world.gc_black:world.gc_dot[1], x-r, y-r, 2*r, 2*r, 0, 360*64);
     }
 
-    if(rend.num == 1){
+    if(rend->num == 1){
       char text[16];
       snprintf(text, sizeof(text), "%d", k+1);
       XDRAWSTRING(world.dis, world.canv, world.gc_black, x, y, text, strlen(text));
     }
-    else if(rend.num == -1){
+    else if(rend->num == -1){
       char text[16];
       const char * s = getname(q);
       s ? snprintf(text, sizeof(text), "%s", s) :  snprintf(text, sizeof(text), "%d", q );
       XDRAWSTRING(world.dis, world.canv, world.gc_black, x, y, text, strlen(text));
     }
 
-    if(rend.bonds>0){
+    if(rend->bonds>0){
       for(int j=k*BONDS_MAX; j<(k+1)*BONDS_MAX; j++){
         int k1 = ac->bonds.a[j];
         if(k1 == -1 ){
@@ -92,7 +92,7 @@ void ac3_draw(atcoord * ac, rendpars rend){
         }
         double dd = BOND_OFFSET * r / sqrt(r2d);
         XDrawLine(world.dis, world.canv, world.gc_black, x+dd*dx, y+dd*dy, x1, y1);
-        if(rend.bonds==2){
+        if(rend->bonds==2){
           char text[16];
           snprintf(text, sizeof(text), "%.3lf", ac->bonds.r[j]);
           XDRAWSTRING(world.dis, world.canv, world.gc_black, x+dx/2, y+dy/2, text, strlen(text));
diff --git a/src/v/ac3_print.c b/src/v/ac3_print.c
index 64a9d1d..b7e22eb 100644
--- a/src/v/ac3_print.c
+++ b/src/v/ac3_print.c
@@ -1,15 +1,16 @@
 #include "v.h"
+#include "matrix.h"
 #include "vec3.h"
 
-void ac3_print(atcoord * ac, rendpars rend){
+void ac3_print(atcoord * ac, rendpars * rend){
   PRINTOUT(stdout, "$molecule\ncart\n");
   for(int k=0; k<ac->n; k++){
     PRINTOUT(stdout, "%3d   % lf   % lf   % lf",
         ac->q[k],
-        rend.xy0[0] + ac->r[k*3  ],
-        rend.xy0[1] + ac->r[k*3+1],
+        rend->xy0[0] + ac->r[k*3  ],
+        rend->xy0[1] + ac->r[k*3+1],
                  ac->r[k*3+2]);
-    if(rend.bonds>0){
+    if(rend->bonds>0){
       for(int j=0; j<BONDS_MAX; j++){
         int k1 = ac->bonds.a[k*BONDS_MAX+j];
         if(k1 == -1 ){
@@ -24,8 +25,20 @@ void ac3_print(atcoord * ac, rendpars rend){
   return;
 }
 
-void ac3_print_xyz(atcoord * ac, rendpars rend){
-  PRINTOUT(stdout, "%d\n\n", ac->n);
+void ac3_print_xyz(atcoord * ac, rendpars * rend){
+  PRINTOUT(stdout, "%d\n", ac->n);
+  if(ac->cell.boundary==CELL){
+    double C[9];
+    mx_multmx(3,3,3, C, rend->ac3rmx, ac->cell.rot_to_lab_basis);
+    PRINTOUT(stdout, "Lattice=\"%lf %lf %lf %lf %lf %lf %lf %lf %lf\"\n",
+        C[0], C[3], C[6],
+        C[1], C[4], C[7],
+        C[2], C[5], C[8]);
+  }
+  else{
+    PRINTOUT(stdout, "\n");
+  }
+
   for(int k=0; k<ac->n; k++){
     const char * s = getname(ac->q[k]);
     int ok = s && s[0];
@@ -36,32 +49,32 @@ void ac3_print_xyz(atcoord * ac, rendpars rend){
       PRINTOUT(stdout, " %3d", ac->q[k]);
     }
     PRINTOUT(stdout, "   % lf   % lf   % lf\n",
-        rend.xy0[0] + ac->r[k*3  ],
-        rend.xy0[1] + ac->r[k*3+1],
-                 ac->r[k*3+2]);
+        rend->xy0[0] + ac->r[k*3  ],
+        rend->xy0[1] + ac->r[k*3+1],
+                       ac->r[k*3+2]);
   }
   return;
 }
 
-void ac3_print2fig(atcoord * ac, rendpars rend, double vert[9]){
+void ac3_print2fig(atcoord * ac, rendpars * rend){
   int n = ac->n;
   for(int i=0; i<n; i++){
     PRINTOUT(stdout, "atom %3d% 13.7lf% 13.7lf% 13.7lf\n", ac->q[i],
-        rend.xy0[0] + ac->r[i*3  ],
-        rend.xy0[1] + ac->r[i*3+1],
+        rend->xy0[0] + ac->r[i*3  ],
+        rend->xy0[1] + ac->r[i*3+1],
                  ac->r[i*3+2]);
   }
 
-  if(vert){
+  if(ac->cell.boundary==CELL){
     for(int i=0; i<8; i++){
       double v[3];
-      r3mx(v, vert+3*i, rend.ac3rmx);
+      r3mx(v, ac->cell.vertices+3*i, rend->ac3rmx);
       PRINTOUT(stdout, "atom %3d% 13.7lf% 13.7lf% 13.7lf\n", 0,
-          rend.xy0[0] + v[0], rend.xy0[1] + v[1], v[2]);
+          rend->xy0[0] + v[0], rend->xy0[1] + v[1], v[2]);
     }
   }
 
-  if(rend.bonds>0){
+  if(rend->bonds>0){
     for(int k=0; k<n; k++){
       for(int j=0; j<BONDS_MAX; j++){
         int k1 = ac->bonds.a[k*BONDS_MAX+j];
@@ -75,7 +88,7 @@ void ac3_print2fig(atcoord * ac, rendpars rend, double vert[9]){
     }
   }
 
-  if(vert){
+  if(ac->cell.boundary==CELL){
 #define LINE(I,J)   PRINTOUT(stdout, "bond %3d %3d % 3d\n", (J)+n+1, (I)+n+1, -1)
     for(int i=0; i<8; i+=2){
       LINE(i,i+1); // || z-axis
diff --git a/src/v/ac3_read.c b/src/v/ac3_read.c
index 18e081b..2e5c1df 100644
--- a/src/v/ac3_read.c
+++ b/src/v/ac3_read.c
@@ -1,8 +1,32 @@
 #include "v.h"
 #include "vecn.h"
 #include "vec3.h"
+#include "matrix.h"
 
-atcoord * atcoord_fill(mol * m0, int b, const geompars geom){
+#define EPS_INV 1e-15
+
+static void cell_fill(cellpars * cell, const double abc[9]){
+  const double * a = abc+0;
+  const double * b = abc+3;
+  const double * c = abc+6;
+  for(int i=0; i<2; i++){
+    for(int j=0; j<2; j++){
+      for(int k=0; k<2; k++){
+        r3sums3(cell->vertices + (i*4+j*2+k)*3, a, i-0.5, b, j-0.5, c, k-0.5);
+      }
+    }
+  }
+  double rot_to_lab_basis[9] = {a[0], b[0], c[0],
+                                a[1], b[1], c[1],
+                                a[2], b[2], c[2]};
+  veccp(9,     cell->rot_to_lab_basis, rot_to_lab_basis);
+  mx_id(3,     cell->rot_to_cell_basis);
+  mx_inv(3, 3, cell->rot_to_cell_basis, rot_to_lab_basis, EPS_INV);
+  cell->boundary = CELL;
+  return;
+}
+
+atcoord * atcoord_fill(mol * m0, const int b, const geompars geom, const double cell[9]){
   int n = m0->n;
 
   size_t q_size = sizeof(int   ) * n;
@@ -47,6 +71,19 @@ atcoord * atcoord_fill(mol * m0, int b, const geompars geom){
     center_mol(n, m->r, geom.center==2 ? m->q : NULL);
   }
   veccp(n*3, m->r0, m->r);
+
+  if(geom.boundary == CELL){
+    cell_fill(&m->cell, geom.cell);
+  }
+  else if(geom.boundary == SHELL){
+    m->cell.vertices[0] =  geom.shell[0];
+    m->cell.vertices[1] =  geom.shell[1];
+    m->cell.boundary = SHELL;
+  }
+  else if(cell && geom.boundary!=CELL_DISABLED){
+    cell_fill(&m->cell, cell);
+  }
+
   return m;
 }
 
@@ -54,10 +91,12 @@ atcoord * ac3_read(readpars read, int b, const geompars geom, format_t * format)
 
   mol * m = NULL;
   FILE * f = read.f;
+  int cell_found = 0;
+  double cell[9] = {};
 
   switch(*format){
     case XYZ:
-      if((m=ac3_read_xyz(f))){
+      if((m=ac3_read_xyz(f, &cell_found, cell))){
         *format = XYZ;
       }
       break;
@@ -72,7 +111,7 @@ atcoord * ac3_read(readpars read, int b, const geompars geom, format_t * format)
       }
       break;
     default:
-      if((m=ac3_read_xyz(f))){
+      if((m=ac3_read_xyz(f, &cell_found, cell))){
         *format = XYZ;
       }
       if(!m){
@@ -92,7 +131,7 @@ atcoord * ac3_read(readpars read, int b, const geompars geom, format_t * format)
   }
   m->name = read.fname;
 
-  atcoord * M = atcoord_fill(m, b, geom);
+  atcoord * M = atcoord_fill(m, b, geom, cell_found ? cell : NULL);
   free(m);
   return M;
 }
diff --git a/src/v/ac3_read_xyz.c b/src/v/ac3_read_xyz.c
index 02a3205..9660352 100644
--- a/src/v/ac3_read_xyz.c
+++ b/src/v/ac3_read_xyz.c
@@ -1,33 +1,85 @@
 #include "v.h"
+#include "vecn.h"
 
-mol * ac3_read_xyz(FILE * f){
+static inline int eat_line(FILE * f){
+  char c;
+  do{
+    c = fgetc(f);
+    if(c==EOF){
+      return -1;
+    }
+  } while(c!='\n');
+  return 0;
+}
 
-  // count atoms
-  int n = 0;
-  if(fscanf(f, "%d", &n) != 1){
-    return NULL;
-  }
+int read_header(FILE * f, double cell[9]){
+
+  const char pattern[] = "Lattice=\"";
 
-  // skip the comment line
-  int c = fgetc(f);
+  int c, i=0, found=0;
   do{
     c = fgetc(f);
     if(c==EOF){
-      GOTOHELL;
+      return -1;
+    }
+    if(c==pattern[i]){
+      if(++i==sizeof(pattern)-1){
+        found = 1;
+        break;
+      }
+    }
+    else{
+      i=0;
     }
   } while(c!='\n');
 
+  if(!found){
+    return 0;
+  }
+
+  int ret = (fscanf(f, "%lf%lf%lf%lf%lf%lf%lf%lf%lf", cell, cell+1, cell+2, cell+3, cell+4, cell+5, cell+6, cell+7, cell+8)==9);
+
+  if(eat_line(f)){
+    return -1;
+  }
+
+  return ret;
+}
+
+mol * ac3_read_xyz(FILE * f, int * cell_found, double * cell){
+
+  // count atoms
+  int n = 0;
+  if(fscanf(f, "%d", &n) != 1){
+    return NULL;
+  }
+  fgetc(f);
+
+  // read header
+  double tcell[9];
+  int r = read_header(f, tcell);
+  if(r==-1){
+    return NULL;
+  }
+  else if(r==1){
+    veccp(9, cell, tcell);
+    *cell_found = 1;
+  }
+
   // fill in
   mol * m = alloc_mol(n);
-  char tmp_str[BIGSTRLEN];
   for(int i=0; i<n; i++){
     styp type;
-    if (fscanf (f, "%7s%lf%lf%lf%[^\n]",
-          type, m->r+3*i, m->r+3*i+1, m->r+3*i+2, tmp_str) < 4) {
+    if (fscanf (f, "%7s%lf%lf%lf", type, m->r+3*i, m->r+3*i+1, m->r+3*i+2) != 4){
       free(m);
       return NULL;
     }
     m->q[i] = get_element(type);
+
+    if(eat_line(f)){
+      free(m);
+      return NULL;
+    }
   }
 
   return m;
diff --git a/src/v/cli.c b/src/v/cli.c
index 21c032c..80a55a9 100644
--- a/src/v/cli.c
+++ b/src/v/cli.c
@@ -1,9 +1,8 @@
 #include "v.h"
 #include "vecn.h"
 #include "matrix.h"
-#include "vec3.h"
 
-#define EPS_INV 1e-15
+#define EPS 1e-15
 
 static int lazysscanf(char * const s, const char * const template, char ** ret){
   // if the string `s` begins with the template,
@@ -37,15 +36,15 @@ static int sscan_rot(const char * arg, double ac3rmx[9]){
   return 1;
 }
 
-static int sscan_cell(const char * arg, cellpars * cp){
+static int sscan_cell(const char * arg, geompars * geom){
   double cell[9];
   int count = sscanf(arg, "cell:b%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf", cell, cell+1, cell+2, cell+3, cell+4, cell+5, cell+6, cell+7, cell+8);
   if(count > 0){
-    if(count==9 || count==3){
+    if(count==9 || count==3 || count==1){
       vecscal(count, cell, BA);
     }
     else{
-      PRINT_WARN("option `cell:b` takes exactly 3 or 9 comma-separated arguments\n")
+      PRINT_WARN("option `cell:b` takes exactly 1, 3 or 9 comma-separated arguments\n")
       return 1;
     }
   }
@@ -54,61 +53,56 @@ static int sscan_cell(const char * arg, cellpars * cp){
     if(count <= 0){
       return 0;
     }
-    else if(count!=3 && count!=9){
-      PRINT_WARN("option `cell:` takes exactly 3 or 9 comma-separated arguments\n")
+    else if(count!=3 && count!=9 && count!=1){
+      PRINT_WARN("option `cell:` takes exactly 1, 3, or 9 comma-separated arguments\n")
       return 1;
     }
   }
 
-  double a[3]={}, b[3]={}, c[3]={};
-  if(count==3){
-    a[0] = cell[0];
-    b[1] = cell[1];
-    c[2] = cell[2];
+  double s=0.0;
+  for(int i=0; i<count; i++){
+    s += fabs(cell[i]);
   }
-  else{
-    r3cp(a, cell+0);
-    r3cp(b, cell+3);
-    r3cp(c, cell+6);
+  if(s<EPS){
+    geom->boundary = CELL_DISABLED;
+    return 1;
   }
 
-  for(int i=0; i<2; i++){
-    for(int j=0; j<2; j++){
-      for(int k=0; k<2; k++){
-        r3sums3(cp->vertices + (i*4+j*2+k)*3, a, i-0.5, b, j-0.5, c, k-0.5);
-      }
-    }
+  if(count==1){
+    geom->cell[0] = cell[0];
+    geom->cell[4] = cell[0];
+    geom->cell[8] = cell[0];
   }
-
-  double rot_to_lab_basis[9] = {a[0], b[0], c[0],
-                                a[1], b[1], c[1],
-                                a[2], b[2], c[2]};
-  veccp(9,     cp->rot_to_lab_basis, rot_to_lab_basis);
-  mx_id(3,     cp->rot_to_cell_basis);
-  mx_inv(3, 3, cp->rot_to_cell_basis, rot_to_lab_basis, EPS_INV);
-
-  cp->vert = 1;
+  else if(count==3){
+    geom->cell[0] = cell[0];
+    geom->cell[4] = cell[1];
+    geom->cell[8] = cell[2];
+  }
+  else{
+    veccp(9, geom->cell, cell);
+  }
+  geom->boundary = CELL;
   return 1;
 }
 
-static int sscan_shell(const char * arg, cellpars * cp){
+static int sscan_shell(const char * arg, geompars * geom){
   double shell[2];
   int count = sscanf(arg, "shell:b%lf,%lf", shell, shell+1);
   if(count > 0){
-    if(count==2)
-      vecscal(count, shell, BA);
+    vecscal(count, shell, BA);
   }
   else{
     count = sscanf(arg, "shell:%lf,%lf", shell, shell+1);
   }
-  if(count <= 0)
+  if(count>0){
+    geom->shell[0] = shell[0];
+    geom->shell[1] = shell[count>1];
+    geom->boundary = SHELL;
+    return 1;
+  }
+  else{
     return 0;
-  if(count == 2){
-    cp->vertices[0] = shell[0];
-    cp->vertices[1] = shell[1];
-    cp->vert = 2;
   }
-  return 1;
 }
 
 static int cli_parse_arg(char * arg, allpars * ap){
@@ -135,8 +129,8 @@ static int cli_parse_arg(char * arg, allpars * ap){
     || lazysscanf(arg, "exitcom:", &dp->ui.on_exit)
     || lazysscanf(arg, "colors:",  &ts)
     || sscan_rot  (arg, dp->rend.ac3rmx)
-    || sscan_cell (arg, &dp->cell)
-    || sscan_shell(arg, &dp->cell)
+    || sscan_cell (arg, &dp->geom)
+    || sscan_shell(arg, &dp->geom)
     ;
 
   if(vib==0){
@@ -194,7 +188,6 @@ static allpars allpars_init(void){
   ap.dp.anim.dt     = DEFAULT_TIMEOUT;
   ap.dp.anal.symtol = DEFAULT_SYMTOL;
   ap.dp.bond.rl     = 1.0;
-  ap.dp.cell.vert   = -1;
   ap.dp.geom.center = 1;
 
   ap.dp.rend.r     = 1.0;
diff --git a/src/v/evr.c b/src/v/evr.c
index f502917..d81ae67 100644
--- a/src/v/evr.c
+++ b/src/v/evr.c
@@ -38,14 +38,21 @@ void kp_readagain(object * ent, drawpars * dp){
 
 void kp_print(object * ent, drawpars * dp){
   if (dp->task == AT3COORDS){
-    ac3_print(ent->m[dp->n], dp->rend);
+    ac3_print(ent->m[dp->n], &dp->rend);
   }
   return;
 }
 
 void kp_print_xyz(object * ent, drawpars * dp){
   if (dp->task == AT3COORDS){
-    ac3_print_xyz(ent->m[dp->n], dp->rend);
+    ac3_print_xyz(ent->m[dp->n], &dp->rend);
+  }
+  return;
+}
+
+void kp_print2fig(object * ent, drawpars * dp){
+  if (dp->task == AT3COORDS){
+    ac3_print2fig(ent->m[dp->n], &dp->rend);
   }
   return;
 }
@@ -60,13 +67,6 @@ void kp_printrot(object * ent __attribute__ ((unused)), drawpars * dp){
   return;
 }
 
-void kp_print2fig(object * ent, drawpars * dp){
-  if (dp->task == AT3COORDS){
-    ac3_print2fig(ent->m[dp->n], dp->rend, dp->cell.vert==1?dp->cell.vertices:NULL);
-  }
-  return;
-}
-
 static void rl_changed(object * ent, drawpars * dp){
   for(int i=0; i<ent->n; i++){
     ent->m[i]->bonds.flag = 0;
@@ -209,38 +209,46 @@ static void mol2cell(double r[3], cellpars * cell){
   return;
 }
 
-static void move_pbc(object * acs, drawpars * dp, int dir, double d){
-
-  double dr[3], v[3] = {};
-  v[dir] = d;  // translation in the view basis
-  r3mxt(dr, v, dp->rend.ac3rmx);  // translation in the mol basis.
-                                  // not true if the initial "rotation" from CLI is not unitary, but ignore this
-  for(int i=0; i<acs->n; i++){
-    atcoord * m = acs->m[i];
-    for(int j=0; j<m->n; j++){
-      double * r = m->r0+j*3;
-      r3add(r, dr);
-      mol2cell(r, &dp->cell);
-      r3cp(m->r+j*3, r);
-      m->rotated = 0;
-    }
-    if(dp->rend.bonds>0){
-      m->bonds.flag = 0;
-      m->bonds.rl *= RL_MOVE_PBC_SCALE;
-    }
+static void move_pbc(atcoord * m, drawpars * dp, double dr[3]){
+  for(int j=0; j<m->n; j++){
+    double * r = m->r0+j*3;
+    r3add(r, dr);
+    mol2cell(r, &m->cell);///TODO
+    r3cp(m->r+j*3, r);
+    m->rotated = 0;
+  }
+  if(dp->rend.bonds>0){
+    m->bonds.flag = 0;
+    m->bonds.rl *= RL_MOVE_PBC_SCALE;
   }
   return;
 }
 
 static void move_ent(object * ent, drawpars * dp, int dir, double step){
+  atcoord * m = ent->m[dp->n];
+
   if(dp->ui.modkey){
     step *= STEP_MOD;
   }
-  if(dp->cell.vert == 1){
-    move_pbc(ent, dp, dir, step);
-  }
-  else {
+
+  if((dp->task==VIBRO) || (m->cell.boundary != CELL)){
     dp->rend.xy0[dir] += step;
+    return;
+  }
+
+  double dr[3], v[3] = {};
+  v[dir] = step;  // translation in the view basis
+  r3mxt(dr, v, dp->rend.ac3rmx);  // translation in the mol basis.
+                                  // not true if the initial "rotation" from CLI is not unitary, but ignore this
+
+  if(dp->geom.boundary==CELL){
+    // all have the same cell -> move together, otherwise move separately
+    for(int i=0; i<ent->n; i++){
+      move_pbc(ent->m[i], dp, dr);
+    }
+  }
+  else{
+    move_pbc(m, dp, dr);
   }
   return;
 }
diff --git a/src/v/load.c b/src/v/load.c
index d211db9..2210c9a 100644
--- a/src/v/load.c
+++ b/src/v/load.c
@@ -191,7 +191,7 @@ object * acs_from_var(int n, mol * m, vibr_t vib, allpars * ap){
     ent->vib = NULL;
 
     for(int i=0; i<n; i++){
-      ent->m[i] = atcoord_fill(m+i, dp->rend.bonds, dp->geom);
+      ent->m[i] = atcoord_fill(m+i, dp->rend.bonds, dp->geom, NULL);
     }
 
     fill_nf(ent, 0);
@@ -207,7 +207,7 @@ object * acs_from_var(int n, mol * m, vibr_t vib, allpars * ap){
   else{
     ent->Nmem = ent->n = 1;
     ent->m     = malloc(ent->Nmem*sizeof(atcoord *));
-    ent->m[0] = atcoord_fill(m+n-1, dp->rend.bonds, dp->geom);
+    ent->m[0] = atcoord_fill(m+n-1, dp->rend.bonds, dp->geom, NULL);
     int nat = ent->m[0]->n;
     ent->vib = make_vibr_t(vib.n, nat);
     veccp(vib.n*nat*3, ent->vib->disp, vib.disp);
diff --git a/src/v/pars.h b/src/v/pars.h
index a536b12..01029d0 100644
--- a/src/v/pars.h
+++ b/src/v/pars.h
@@ -9,6 +9,13 @@ typedef enum {
   CPK_COLORS,
 } colorscheme_t;
 
+typedef enum {
+  NO_BOUNDARY,
+  CELL,
+  SHELL,
+  CELL_DISABLED,
+} boundary_t;
+
 typedef struct {
   colorscheme_t  colors;   // colorscheme (v or cpk)
   int               gui;   // if gui is enabled
@@ -26,6 +33,10 @@ typedef struct {
   int center;   // 0: nothing; 1: center each molecule upon reading; 2: center wrt center of mass
   int inertia;  // 0: nothing; 1: rotate each molecule upon reading wrt axis of inertia
   int bohr;     // 0: Å        1: Bohr
+
+  boundary_t boundary;
+  double cell[9];
+  double shell[2];
 } geompars;
 
 typedef struct {
@@ -55,7 +66,7 @@ typedef struct {
   double vertices[3*8];           // parameters of cell/shell
   double rot_to_lab_basis[3*3];   // "rotation" matrix for PBC
   double rot_to_cell_basis[3*3];  // "rotation" matrix for PBC
-  int    vert;          // 0: nothing;     1: show cell;      2: show shell
+  boundary_t boundary;
 } cellpars;
 
 typedef struct {
@@ -81,7 +92,6 @@ typedef struct {
   uipars    ui;
   rendpars  rend;
   bondpars  bond;
-  cellpars  cell;
   animpars  anim;
   analpars  anal;
 
diff --git a/src/v/redraw.c b/src/v/redraw.c
index 190b134..a077034 100644
--- a/src/v/redraw.c
+++ b/src/v/redraw.c
@@ -1,6 +1,18 @@
 #include "v.h"
 #include "evr.h"
 
+static void draw_boundary(atcoord * m , rendpars * rend){
+  if(m->cell.boundary==CELL){
+    double v[8*3];
+    rot3d(8, v, m->cell.vertices, rend->ac3rmx);
+    draw_vertices(v, rend);
+  }
+  else if(m->cell.boundary==SHELL){
+    draw_shell(m->cell.vertices, rend);
+  }
+  return;
+}
+
 static void ac3_text(atcoord * ac, drawpars * dp){
   char text[32];
   char text_fname[STRLEN];
@@ -79,34 +91,30 @@ void redraw_ac3(object * ent, drawpars * dp){
   atcoord * m = ent->m[dp->n];
   fill_bonds(m, dp);
   rotate_mol(m, dp);
+  if(m->cell.boundary==CELL){
+    dp->rend.xy0[0] = 0.0;
+    dp->rend.xy0[1] = 0.0;
+  }
 
   clear_canv();
-  ac3_draw(m, dp->rend);
+  ac3_draw(m, &dp->rend);
   ac3_text(m, dp);
-  if(dp->cell.vert == 1){
-    double v[8*3];
-    rot3d(8, v, dp->cell.vertices, dp->rend.ac3rmx);
-    drawvertices(v, dp->rend);
-  }
-  else if(dp->cell.vert == 2){
-    drawshell(dp->cell.vertices, dp->rend);
-  }
+  draw_boundary(m, &dp->rend);
   fill_canv();
   return;
 }
 
 void redraw_vibro(object * ent, drawpars * dp){
-
   atcoord * m  = ent->m[0];
-  double  * dr = ent->vib->disp + dp->n * m->n*3;
-
   fill_bonds(m, dp);
+  double  * dr = ent->vib->disp + dp->n * m->n*3;
   vecsums(m->n*3, m->r, m->r0, dr, VIBR_AMP*sqrt(m->n)*sin(dp->anim.t * 2.0*M_PI/TMAX));
   rot3d_inplace(m->n, m->r, dp->rend.ac3rmx);
 
   clear_canv();
-  ac3_draw(m, dp->rend);
+  ac3_draw(m, &dp->rend);
   vibro_text(ent->vib, dp);
+  draw_boundary(m, &dp->rend);
   fill_canv();
   return;
 }
diff --git a/src/v/scale.c b/src/v/scale.c
index 4721933..a5ba651 100644
--- a/src/v/scale.c
+++ b/src/v/scale.c
@@ -15,6 +15,17 @@ double ac3_scale(atcoord * ac){
     double d2  = r3d2(center, ac->r+3*k);
     d2max = MAX(d2max, d2+rad*rad);
   }
+  if(ac->cell.boundary==SHELL){
+    for(int i=0; i<2; i++){
+      d2max = MAX(d2max, 2*ac->cell.vertices[i]*ac->cell.vertices[i]);
+    }
+  }
+  else if(ac->cell.boundary==CELL){
+    for(int k=0; k<8; k++){
+      double d2 = r3d2(center, ac->cell.vertices+3*k);
+      d2max = MAX(d2max, d2);
+    }
+  }
   return VIEWPORT_FILL / sqrt(d2max);
 }
 
diff --git a/src/v/v.h b/src/v/v.h
index 05aecae..f8dc1ff 100644
--- a/src/v/v.h
+++ b/src/v/v.h
@@ -6,7 +6,6 @@
 #define BONDS_MAX 32
 #define POINTER_SPEED 2.0
 #define STRLEN 256
-#define BIGSTRLEN 4096
 #define MAX_LINES 5
 
 #include "pars.h"
@@ -39,6 +38,7 @@ typedef struct {
   int nf[2];             // number of molecule in file, file size
   styp   sym;            // point group
   bondstr bonds;
+  cellpars cell;
 } atcoord;
 
 typedef struct {
@@ -68,11 +68,11 @@ vibr_t * make_vibr_t(int n_modes, int n_atoms);
 vibr_t * mode_read(FILE * f, int na);
 // ac3_read*.c
 int read_cart_atom(FILE * f, int n, mol * m);
-atcoord * atcoord_fill(mol * m, int b, geompars geom);
+atcoord * atcoord_fill(mol * m0, const int b, const geompars geom, const double cell[9]);
 atcoord * ac3_read(readpars read, int b, geompars geom, format_t * format);
 mol * ac3_read_in (FILE * f);
 mol * ac3_read_out(FILE * f);
-mol * ac3_read_xyz(FILE * f);
+mol * ac3_read_xyz(FILE * f, int * lattice_found, double * lattice);
 
 // man.c
 void printman(FILE * f, char * exename);
@@ -87,11 +87,11 @@ void redraw_ac3(object * ent, drawpars * dp);
 void redraw_vibro(object * ent, drawpars * dp);
 
 // ac3_draw.c
-void ac3_draw      (atcoord * ac, rendpars rend);
+void ac3_draw     (atcoord * ac, rendpars * rend);
 // ac3_print.c
-void ac3_print    (atcoord * ac, rendpars rend);
-void ac3_print_xyz(atcoord * ac, rendpars rend);
-void ac3_print2fig(atcoord * ac, rendpars rend, double * v);
+void ac3_print    (atcoord * ac, rendpars * rend);
+void ac3_print_xyz(atcoord * ac, rendpars * rend);
+void ac3_print2fig(atcoord * ac, rendpars * rend);
 // bonds.c
 void bonds_fill(bondpars bond, atcoord * ac);
 
@@ -107,8 +107,8 @@ void init_x       (const char * const capt, const colorscheme_t colorscheme);
 void init_font    (char * fontname);
 void textincorner (const char * const lines[MAX_LINES], const int red[MAX_LINES]);
 void setcaption   (const char * const capt);
-void drawvertices (double * v, rendpars rend);
-void drawshell    (double r[2], rendpars rend);
+void draw_vertices(double * v, rendpars * rend);
+void draw_shell   (double r[2], rendpars * rend);
 int  savepic      (char * s);
 void clear_canv();
 void fill_canv();
diff --git a/src/v/x.c b/src/v/x.c
index d28333a..b6df6e2 100644
--- a/src/v/x.c
+++ b/src/v/x.c
@@ -153,14 +153,14 @@ void setcaption(const char * const capt){
   return;
 }
 
-void draw_edge(double vi[3], double vj[3], rendpars rend){
+void draw_edge(double vi[3], double vj[3], rendpars * rend){
   int iw = (vi[2]>0.0 || vj[2]>0.0) ? 0 : 1;
   XDrawLine(world.dis, world.canv, world.gc_dot[iw],
       SCREEN_X(vi[0]), SCREEN_Y(vi[1]), SCREEN_X(vj[0]), SCREEN_Y(vj[1]));
   return;
 }
 
-void drawvertices(double * v, rendpars rend){
+void draw_vertices(double * v, rendpars * rend){
 #define LINE(i,j) draw_edge(v+(i)*3, v+(j)*3, rend)
   for(int i=0; i<8; i+=2){
     LINE(i,i+1); // || z-axis
@@ -175,8 +175,8 @@ void drawvertices(double * v, rendpars rend){
   return;
 }
 
-void drawshell(double r[2], rendpars rend){
-  double d = world.size * rend.scale;
+void draw_shell(double r[2], rendpars * rend){
+  double d = world.size * rend->scale;
   for(int i=0; i<2; i++){
     XDrawArc(world.dis, world.canv, world.gc_dot[1-i],
         SCREEN_X(-r[i]), SCREEN_Y(r[i]),
diff --git a/src/v/x.h b/src/v/x.h
index 17698c2..e85aa30 100644
--- a/src/v/x.h
+++ b/src/v/x.h
@@ -7,8 +7,8 @@
 #define NCOLORS 110
 #define LINE_WIDTH 2
 
-#define SCREEN_X(X)  (world.W/2 + world.size * rend.scale*(rend.xy0[0] + (X)))
-#define SCREEN_Y(Y)  (world.H/2 - world.size * rend.scale*(rend.xy0[1] + (Y)))
+#define SCREEN_X(X)  (world.W/2 + world.size * rend->scale*(rend->xy0[0] + (X)))
+#define SCREEN_Y(Y)  (world.H/2 - world.size * rend->scale*(rend->xy0[1] + (Y)))
 
 typedef struct {
   Display * dis;