diff --git a/.gitignore b/.gitignore
index c07309633e9f7eb78b1f60a5ce45a617a8a79f8c..a24acea524b5f4ab694b591696d46cf7c6761af2 100644
--- a/.gitignore
+++ b/.gitignore
@@ -19,6 +19,7 @@
 # ViM, emacs, nano, leafpad
 *.swp
 *~
+tags
 
 # Python
 *.pyc
@@ -54,7 +55,12 @@ stamp-h1
 *.tab.c
 *.tab.h
 
+# PBS/torque job logs
+*.o[0-9]*
+*.e[0-9]*
+
 # General
 *.tar.gz
 *.pdf
 */bin
+*.txt
diff --git a/data/benchmarks/water/ffield_Achtyl b/data/benchmarks/water/ffield_Achtyl
new file mode 100644
index 0000000000000000000000000000000000000000..ce54cb3bad2a701ff11a41f08679d98026a022d9
--- /dev/null
+++ b/data/benchmarks/water/ffield_Achtyl
@@ -0,0 +1,424 @@
+Reactive MD-force field: water/C/Ca/Ti AKS2 July 18; water from Feb 29          
+ 39       ! Number of general parameters                                        
+   50.0000 !Overcoordination parameter                                          
+    9.5469 !Overcoordination parameter                                          
+    1.6725 !Valency angle conjugation parameter                                 
+    1.7224 !Triple bond stabilisation parameter                                 
+    6.8702 !Triple bond stabilisation parameter                                 
+   54.6742 !C2-correction                                                       
+    1.0588 !Undercoordination parameter                                         
+    4.6000 !Triple bond stabilisation parameter                                 
+   12.1176 !Undercoordination parameter                                         
+   13.3056 !Undercoordination parameter                                         
+  -34.5448 !Triple bond stabilization energy                                    
+    0.0000 !Lower Taper-radius                                                  
+   10.0000 !Upper Taper-radius                                                  
+    2.8793 !Not used                                                            
+   33.8667 !Valency undercoordination                                           
+    6.0891 !Valency angle/lone pair parameter                                   
+    1.0563 !Valency angle                                                       
+    2.0384 !Valency angle parameter                                             
+    6.1431 !Not used                                                            
+    6.9290 !Double bond/angle parameter                                         
+    0.3989 !Double bond/angle parameter: overcoord                              
+    3.9954 !Double bond/angle parameter: overcoord                              
+   -2.4837 !Not used                                                            
+    8.5385 !Torsion/BO parameter                                                
+    6.7491 !Torsion overcoordination                                            
+    0.1414 !Torsion overcoordination                                            
+   -1.2327 !Conjugation 0 (not used)                                            
+    1.1348 !Conjugation                                                         
+    1.5591 !vdWaals shielding                                                   
+    0.1000 !Cutoff for bond order (*100)                                        
+    1.7602 !Valency angle conjugation parameter                                 
+    0.6991 !Overcoordination parameter                                          
+   50.0000 !Overcoordination parameter                                          
+    1.8512 !Valency/lone pair parameter                                         
+  548.6451 !Softness                                                            
+    0.0000 !Cutoff                                                              
+    5.0000 !Molecular energy (not used)                                         
+    0.0000 !Molecular energy (not used)                                         
+    0.7903 !Valency angle conjugation parameter                                 
+ 13   ! Nr of atoms; cov.r; valency;a.m;Rvdw;Evdw;gammaEEM;cov.r2;#             
+            alfa;gammavdW;valency;Eunder;Eover;chiEEM;etaEEM;n.u.               
+            cov r3;Elp;Heat inc.;n.u.;n.u.;n.u.;n.u.                            
+            ov/un;val1;n.u.;val3,vval4                                          
+ C    1.3417   4.0000  12.0000   2.0464   0.1064   0.4670   1.1681   4.0000     
+      9.0000   1.5000   4.0000  40.0000  79.5548   5.3422   4.5000   0.0000     
+      1.1437   0.0000 181.0000   5.4236  19.2788  13.5366   3.1838   0.0000     
+     -4.0998   4.8750   1.0564   4.0000   2.9663   1.2500   0.4000  21.3612     
+ H    0.9102   1.0000   1.0080   1.0996   0.1000   0.6683  -0.1000   1.0000     
+      9.2193   5.5154   1.0000   0.0000 121.1250   4.9673   6.2079   0.0000     
+     -0.1000   0.0000  55.5000   1.1004   6.8959   0.0003   3.4114   0.0000     
+     -6.5532   3.5000   1.0338   1.0000   2.8793   0.5000   0.1000   1.0000     
+ O    1.2218   2.0000  15.9990   2.2033   0.2841   0.9750   1.0727   6.0000     
+      8.8250 200.0000   4.0000  37.5000 116.0768   8.5000   7.9071   0.0000     
+      0.9049   0.5488  68.0152   2.1943   2.3055   0.0021   5.4479   0.0000     
+     -6.0011   3.6068   1.0493   4.0000   2.9225   0.5000   0.1000   1.0000     
+ N    1.2333   3.0000  14.0000   2.2375   0.1447   1.0000   1.1748   5.0000     
+      9.8626  12.6599   4.0000  30.3181 100.0000   6.0111   6.7037   0.0000     
+      1.0433   0.7872 119.9837   0.7425   6.7920   2.7271   2.2882   0.0000     
+     -2.0000   4.0000   1.0183   4.0000   2.8793   1.0000   0.1000   1.0000     
+ S    1.9405   2.0000  32.0600   2.0677   0.2099   1.0336   1.5479   6.0000     
+      9.9575   4.9055   4.0000  52.9998 112.1416   6.5000   8.2545   0.0000     
+      1.4601   9.7177  71.1843   5.7487  23.2859  12.7147   2.2882   0.0000     
+    -11.0000   2.7466   1.0338   6.2998   2.8793   1.0000   0.1000  10.0000     
+ Mg   1.8315   2.0000  24.3050   2.2464   0.1806   0.5020   1.0000   2.0000     
+     10.9186  27.1205   3.0000  38.0000   0.0000   0.9499   5.6130   0.0000     
+     -1.3000   0.0000 127.9160  49.9248   0.3370   0.0000   2.2882   0.0000     
+     -1.0823   2.3663   1.0564   6.0000   2.9663   1.0000   0.1000  10.0000     
+ P    1.5994   3.0000  30.9738   1.7000   0.1743   1.0000   1.3000   5.0000     
+      9.1909  14.9482   5.0000   0.0000   0.0000   1.6676   7.0946   0.0000     
+     -1.0000  25.0000 125.6300   0.2187  21.4305  15.1425   2.2882   0.0000     
+     -3.9294   3.4831   1.0338   5.0000   2.8793   1.0000   0.1000  10.0000     
+ Na   1.8000   1.0000  22.9898   2.8270   0.1872   0.4000  -1.0000   1.0000     
+     10.0000   2.5000   1.0000   0.0000   0.0000  -0.9871   6.7728   0.0000     
+     -1.0000   0.0000  23.0445 100.0000   1.0000   0.0000   2.2882   0.0000     
+     -2.5000   3.9900   1.0338   8.0000   2.5791   1.0000   0.1000  10.0000     
+ Ti   2.0254   4.0000  47.8800   2.2105   0.1574   0.6311   0.1000   4.0000     
+     12.7041  16.6482   4.0000   0.1000   0.0000  -1.3647   6.8406   0.0000     
+     -1.0000   0.0000 143.1770  27.6505  -0.0753   0.0064   2.2882   0.0000     
+    -15.0000   3.8359   1.0338  12.0000   2.2632   1.0000   0.1000  10.0000     
+ Cl   1.7140   1.0000  35.4500   1.9139   0.2000   0.3500  -1.0000   7.0000     
+     11.5345  10.1330   1.0000   0.0000   0.0000   9.9704   6.1703   0.0000     
+     -1.0000   1.2769 143.1770   6.2293   5.2294   0.1542   2.2882   0.0000     
+    -10.2080   2.9867   1.0338   6.2998   2.5791   1.0000   0.1000  10.0000     
+ F    1.2100   1.0000  18.9984   1.8601   0.1200   0.3000  -0.1000   7.0000     
+     11.5000   7.5000   4.0000   9.2533   0.2000   9.0000  15.0000   0.0000     
+     -1.0000  35.0000   1.5000   6.9821   4.1799   1.0561   2.2882   0.0000     
+     -7.3000   2.6656   1.0493   4.0000   2.9225   1.0000   0.1000  10.0000     
+ Ca   1.9058   2.0000  40.0870   2.3698   0.3719   0.6038  -1.0000   2.0000     
+      9.8681   5.0000   3.0000  38.0000   0.0000  -5.7471   7.3556   0.0000     
+     -1.3000   0.0000 220.0000  49.9248   0.3370   0.0000   3.7000   0.0000     
+     -3.9745   3.0069   1.0564   8.0000   2.9663   0.5000   0.1000  10.0000     
+ X   -0.1000   2.0000   1.0080   2.0000   0.0000   1.0000  -0.1000   6.0000     
+     10.0000   2.5000   4.0000   0.0000   0.0000  -0.1000  25.0000   0.0000     
+     -0.1000   0.0000  -2.3700   8.7410  13.3640   0.6690   0.1000   0.0000     
+    -11.0000   2.7466   1.0338   4.0000   2.8793   1.0000   0.1000  10.0000     
+ 62     ! Nr of bonds; Edis1;LPpen;n.u.;pbe1;pbo5;13corr;pbo6                   
+                         pbe2;pbo3;pbo4;Etrip;pbo1;pbo2;ovcorr                  
+  1  1  78.4266 115.3834  68.1631   0.5777  -0.2901   1.0000  34.9989   0.5101  
+         3.8560  -0.1640   8.2326   1.0000  -0.0585   6.7997   1.0000   0.0000  
+  1  2 191.4398   0.0000   0.0000  -0.6539   0.0000   1.0000   6.0000   0.4898  
+         6.3962   1.0000   0.0000   1.0000  -0.0607   6.9960   0.0000   0.0000  
+  2  2 161.2898   0.0000   0.0000  -0.2387   0.0000   1.0000   6.0000   0.6279  
+        13.2089   1.0000   0.0000   1.0000  -0.1674   6.9118   0.0000   0.0000  
+  1  3 170.8787 152.2997 125.6008   0.2914  -0.2305   1.0000  16.7601   0.8645  
+         3.9615  -0.5703   7.4065   1.0000  -0.2514   4.5085   0.0000   0.0000  
+  3  3 222.6268  93.7249  50.8293   0.5588  -0.1000   1.0000  29.7503   0.0125  
+         0.5865  -0.2030   9.1777   1.0000  -0.1357   6.4747   1.0000   0.0000  
+  1  4 163.8300 145.4458  89.6879  -1.3368  -0.3468   1.0000  27.5160   0.1575  
+         0.1817  -0.3114   7.1789   1.0000  -0.2345   4.5111   1.0000   0.0000  
+  3  4 130.8596 169.4551  40.0000   0.3837  -0.1639   1.0000  35.0000   0.2000  
+         1.0000  -0.3579   7.0004   1.0000  -0.1193   6.8773   1.0000   0.0000  
+  4  4 157.9384  82.5526 152.5336   0.4010  -0.1034   1.0000  12.4261   0.5828  
+         0.1578  -0.1509  11.9186   1.0000  -0.0861   5.4271   1.0000   0.0000  
+  2  3 143.3409   0.0000   0.0000  -0.4886   0.0000   1.0000   6.0000   0.3927  
+         1.3588   1.0000   0.0000   0.0000  -0.0781   4.0330   0.0000   0.0000  
+  2  4 210.1187   0.0000   0.0000  -0.3705   0.0000   1.0000   6.0000   0.3284  
+         5.8196   1.0000   0.0000   1.0000  -0.1104   5.5184   0.0000   0.0000  
+  1  5 128.9942  74.5848  55.2528   0.1035  -0.5211   1.0000  18.9617   0.6000  
+         0.2949  -0.2398   8.1175   1.0000  -0.1029   5.6731   1.0000   0.0000  
+  2  5 151.5159   0.0000   0.0000  -0.4721   0.0000   1.0000   6.0000   0.6000  
+         9.4366   1.0000   0.0000   1.0000  -0.0290   7.0050   1.0000   0.0000  
+  3  5   0.0000   0.0000   0.0000   0.5563  -0.4038   1.0000  49.5611   0.6000  
+         0.4259  -0.4577  12.7569   1.0000  -0.1100   7.1145   1.0000   0.0000  
+  4  5   0.0000   0.0000   0.0000   0.4438  -0.2034   1.0000  40.3399   0.6000  
+         0.3296  -0.3153   9.1227   1.0000  -0.1805   5.6864   1.0000   0.0000  
+  5  5  96.1871  93.7006  68.6860   0.0955  -0.4781   1.0000  17.8574   0.6000  
+         0.2723  -0.2373   9.7875   1.0000  -0.0950   6.4757   1.0000   0.0000  
+  2  6  58.6896   0.0000   0.0000  -0.0203  -0.1418   1.0000  13.1260   0.0230  
+         8.2136  -0.1310   0.0000   1.0000  -0.2692   6.4254   0.0000  24.4461  
+  3  6  87.0227   0.0000  43.3991   0.0030  -0.3000   1.0000  36.0000   0.0250  
+         0.0087  -0.2500  12.0000   1.0000  -0.0439   6.6073   1.0000  24.4461  
+  6  6  32.3808   0.0000   0.0000  -0.0076  -0.2000   0.0000  16.0000   0.2641  
+         4.8726  -0.2000  10.0000   1.0000  -0.0729   4.6319   0.0000   0.0000  
+  1  7 110.0000  92.0000   0.0000   0.2171  -0.1418   1.0000  13.1260   0.6000  
+         0.3601  -0.1310  10.7257   1.0000  -0.0869   5.3302   1.0000   0.0000  
+  2  7   0.1466   0.0000   0.0000   0.2250  -0.1418   1.0000  13.1260   0.6000  
+         0.3912  -0.1310   0.0000   1.0000  -0.1029   9.3302   0.0000   0.0000  
+  3  7 201.0058 194.1410   0.0000   1.0000  -0.5000   1.0000  25.0000   0.4873  
+         0.4358  -0.1571  15.8745   1.0000  -0.2431   6.3823   1.0000   0.0000  
+  4  7 130.0000   0.0000   0.0000   0.2171  -0.1418   1.0000  13.1260   0.6000  
+         0.3601  -0.1310  10.7257   1.0000  -0.0869   5.3302   1.0000   0.0000  
+  6  7   0.1000   0.0000   0.0000   0.2500  -0.5000   1.0000  35.0000   0.6000  
+         0.5000  -0.5000  20.0000   1.0000  -0.2000  10.0000   1.0000   0.0000  
+  7  7   0.0000   0.0000   0.0000   0.2171  -0.5000   1.0000  35.0000   0.6000  
+         0.5000  -0.5000  20.0000   1.0000  -0.2000  10.0000   1.0000   0.0000  
+  1  8   0.0000   0.0000   0.0000  -1.0000  -0.3000   1.0000  36.0000   0.7000  
+        10.1151  -0.3500  25.0000   1.0000  -0.1053   8.2003   1.0000   0.0000  
+  2  8   0.0000   0.0000   0.0000  -1.0000  -0.3000   1.0000  36.0000   0.7000  
+        10.1151  -0.3500  25.0000   1.0000  -0.1053   8.2003   1.0000   0.0000  
+  3  8  45.8933   0.0000   0.0000  -0.1511  -0.3000   1.0000  36.0000   0.3105  
+         5.8448  -0.3500  25.0000   1.0000  -0.0659   7.9140   1.0000   0.0000  
+  4  8   0.0000   0.0000   0.0000  -1.0000  -0.3000   1.0000  36.0000   0.7000  
+        10.1151  -0.3500  25.0000   1.0000  -0.1053   8.2003   1.0000   0.0000  
+  5  8   0.0000   0.0000   0.0000  -1.0000  -0.3000   1.0000  36.0000   0.7000  
+        10.1151  -0.3500  25.0000   1.0000  -0.1053   8.2003   1.0000   0.0000  
+  6  8   0.0000   0.0000   0.0000   0.2500  -0.5000   1.0000  35.0000   0.6000  
+         0.5000  -0.5000  20.0000   1.0000  -0.2000  10.0000   1.0000   0.0000  
+  7  8   0.0000   0.0000   0.0000   0.2500  -0.5000   1.0000  35.0000   0.6000  
+         0.5000  -0.5000  20.0000   1.0000  -0.2000  10.0000   1.0000   0.0000  
+  8  8  64.4508   0.0000   0.0000  -0.3738   0.3000   0.0000  25.0000   0.2158  
+         0.9915  -0.4000  12.0000   1.0000  -0.0515   5.0000   0.0000   0.0000  
+  4  6  50.0000  10.0901   0.0000  -1.0000  -0.3000   1.0000  36.0000   0.7058  
+         0.8567  -0.3487  17.4990   1.0000  -0.0794   8.2232   1.0000   0.0000  
+  1  9   0.0000   0.0000   0.0000  -0.2872  -0.3000   1.0000  36.0000   0.0082  
+         1.7973  -0.2500  20.0000   1.0000  -0.2578   6.5219   1.0000   0.0000  
+  2  9   0.0000   0.0000   0.0000  -0.2872  -0.3000   1.0000  36.0000   0.0082  
+         1.7973  -0.2500  20.0000   1.0000  -0.2578   6.5219   1.0000   0.0000  
+  3  9 130.5629  37.6984   0.0000   0.9228  -0.3000   0.0000  36.0000   0.0850  
+         0.1150  -0.2818  16.1571   1.0000  -0.1343   6.8264   0.0000   0.0000  
+  4  9 130.5629  37.6984   0.0000   0.9228  -0.3000   0.0000  36.0000   0.0850  
+         0.1150  -0.2818  16.1571   1.0000  -0.1343   6.8264   0.0000   0.0000  
+  5  9   0.0000   0.0000   0.0000  -0.2872  -0.3000   1.0000  36.0000   0.0082  
+         1.7973  -0.2500  20.0000   1.0000  -0.2578   6.5219   1.0000   0.0000  
+  6  9   0.0000   0.0000   0.0000  -0.2872  -0.3000   1.0000  36.0000   0.0082  
+         1.7973  -0.2500  20.0000   1.0000  -0.2578   6.5219   1.0000   0.0000  
+  7  9   0.0000   0.0000   0.0000  -0.2872  -0.3000   1.0000  36.0000   0.0082  
+         1.7973  -0.2500  20.0000   1.0000  -0.2578   6.5219   1.0000   0.0000  
+  8  9   0.1000   0.0000   0.0000   0.2500  -0.5000   1.0000  35.0000   0.6000  
+         0.5000  -0.5000  20.0000   1.0000  -0.2000  10.0000   1.0000   0.0000  
+  9  9  80.1930   0.0000   0.0000  -0.8469  -0.2000   0.0000  16.0000   0.2022  
+         0.7528  -0.1924  14.9725   1.0000  -0.0885   5.0000   0.0000   0.0000  
+  1 10   0.0000   0.0000   0.0000   0.5000  -0.2000   0.0000  16.0000   0.5000  
+         1.0001  -0.2000  15.0000   1.0000  -0.1000  10.0000   0.0000   0.0000  
+  2 10  98.9788   0.0000   0.0000  -0.0572  -0.2000   0.0000  16.0000   1.1523  
+         2.2822  -0.2000  15.0000   1.0000  -0.1093   5.1686   0.0000   0.0000  
+  3 10   0.0000   0.0000   0.0000   0.5000  -0.2000   0.0000  16.0000   0.5000  
+         1.0001  -0.2000  15.0000   1.0000  -0.1000  10.0000   0.0000   0.0000  
+  4 10   0.0000   0.0000   0.0000   0.5000  -0.2000   0.0000  16.0000   0.5000  
+         1.0001  -0.2000  15.0000   1.0000  -0.1000  10.0000   0.0000   0.0000  
+  5 10   0.0000   0.0000   0.0000   0.5000  -0.2000   0.0000  16.0000   0.5000  
+         1.0001  -0.2000  15.0000   1.0000  -0.1000  10.0000   0.0000   0.0000  
+  6 10   0.0000   0.0000   0.0000   0.5000  -0.2000   0.0000  16.0000   0.5000  
+         1.0001  -0.2000  15.0000   1.0000  -0.1000  10.0000   0.0000   0.0000  
+  7 10   0.0000   0.0000   0.0000   0.5000  -0.2000   0.0000  16.0000   0.5000  
+         1.0001  -0.2000  15.0000   1.0000  -0.1000  10.0000   0.0000   0.0000  
+  8 10   0.0000   0.0000   0.0000   0.5000  -0.2000   0.0000  16.0000   0.5000  
+         1.0001  -0.2000  15.0000   1.0000  -0.1000  10.0000   0.0000   0.0000  
+  9 10   0.0000   0.0000   0.0000   0.5000  -0.2000   0.0000  16.0000   0.5000  
+         1.0001  -0.2000  15.0000   1.0000  -0.1000  10.0000   0.0000   0.0000  
+ 10 10   0.2500   0.0000   0.0000   0.1803  -0.2000   0.0000  16.0000   0.3356  
+         0.9228  -0.2000  15.0000   1.0000  -0.1178   5.6715   0.0000   0.0000  
+  1 11 237.8781   0.0000   0.0000  -0.7438  -0.5000   1.0000  35.0000   1.0460  
+         3.6661  -0.2500  15.0000   1.0000  -0.0800   5.4719   1.0000   0.0000  
+  2 11   0.0000   0.0000   0.0000  -0.4643   0.0000   1.0000   6.0000   0.6151  
+        12.3710   1.0000   0.0000   1.0000  -0.1008   8.5980   0.0000   0.0000  
+  3 11   0.0000   0.0000   0.0000  -0.4643   0.0000   1.0000   6.0000   0.6151  
+        12.3710   1.0000   0.0000   1.0000  -0.1008   8.5980   0.0000   0.0000  
+  4 11   0.0000   0.0000   0.0000  -0.4643   0.0000   1.0000   6.0000   0.6151  
+        12.3710   1.0000   0.0000   1.0000  -0.1008   8.5980   0.0000   0.0000  
+  5 11   0.0000   0.0000   0.0000  -0.4643   0.0000   1.0000   6.0000   0.6151  
+        12.3710   1.0000   0.0000   1.0000  -0.1008   8.5980   0.0000   0.0000  
+ 11 11 250.0765   0.0000   0.0000   0.2298  -0.3500   1.0000  25.0000   0.8427  
+         0.1167  -0.2500  15.0000   1.0000  -0.1506   7.3516   1.0000   0.0000  
+  1 12   0.0000   0.0000   0.0000  -0.0203  -0.1418   1.0000  13.1260   0.0230  
+         8.2136  -0.2500  20.0000   1.0000  -0.2692   6.4254   0.0000  24.4461  
+  2 12   0.0000   0.0000   0.0000  -0.0203  -0.1418   1.0000  13.1260   0.0230  
+         8.2136  -0.2500  20.0000   1.0000  -0.2692   6.4254   0.0000  24.4461  
+  3 12  49.4055   0.0000   0.0000   0.9603  -0.3000   0.0000  36.0000   0.0025  
+         0.4232  -0.2500  12.0000   1.0000  -0.1619   9.6512   0.0000  24.4461  
+ 12 12  22.8272   0.0000   0.0000   0.6166  -0.2000   0.0000  16.0000   0.8225  
+         1.0000  -0.2000  10.0000   1.0000  -0.0831   4.2291   0.0000   0.0000  
+ 26    ! Nr of off-diagonal terms; Ediss;Ro;gamma;rsigma;rpi;rpi2               
+  1  2   0.1240   1.6326   9.8721   1.1578  -1.0000  -1.0000                    
+  2  3   0.0295   1.3181  10.1225   0.9069  -1.0000  -1.0000                    
+  2  4   0.1294   1.3025   9.8751   1.0415  -1.0000  -1.0000                    
+  1  3   0.2287   2.0452   9.0166   1.4046   1.1866   1.0375                    
+  1  4   0.2000   1.8828   9.7673   1.3387   1.2578   1.1539                    
+  3  4   0.1001   2.3274   9.0974   1.5236   1.0493   1.2531                    
+  1  5   0.1408   1.8161   9.9393   1.7986   1.3021   1.4031                    
+  2  5   0.0895   1.6239  10.0104   1.4640  -1.0000  -1.0000                    
+  3  5   0.1022   1.9887  10.0605   1.5799   1.4000  -1.0000                    
+  4  5   0.1505   1.9000  10.5104   1.8000   1.4000  -1.0000                    
+  2  6   0.0100   1.6000  13.2979   1.8670  -1.0000  -1.0000                    
+  3  6   0.0809   1.7000  11.4606   1.5177  -1.0000  -1.0000                    
+  3  7   0.0534   1.7520  10.4281   1.8000   1.4498  -1.0000                    
+  6  7   0.1801   1.8566   9.8498   0.1000  -1.0000  -1.0000                    
+  3  8   0.0825   1.5904  11.3396   1.5905  -1.0000  -1.0000                    
+  2  9   0.1750   1.7939  13.5000   0.0100  -1.0000  -1.0000                    
+  3  9   0.1200   1.8000  10.5000   1.6526   1.4718  -1.0000                    
+  1  9   0.2950   2.2000  11.0937   0.0100  -1.0000  -1.0000                    
+  2 10   0.0376   1.6671   9.6285   1.2123  -1.0000  -1.0000                    
+  3 10   0.1945   2.2766  11.2353  -1.0000  -1.0000  -1.0000                    
+  1 11   0.1071   1.6243  11.0402   1.3176  -1.0000  -1.0000                    
+  2 11   0.0431   1.7204  10.3632   0.5386  -1.0000  -1.0000                    
+  4  9   0.1200   1.8000  10.5000   1.6526   1.4718  -1.0000                    
+  1 12   0.2000   1.5000  14.0000   0.0010   0.0010  -1.0000                    
+  2 12   0.0100   1.0610   9.7343   0.0010   0.0010  -1.0000                    
+  3 12   0.1515   1.8913  12.5160   2.0022  -1.0000  -1.0000                    
+ 98    ! Nr of angles;at1;at2;at3;Thetao,o;ka;kb;pv1;pv2                        
+  1  1  1  74.4118  38.0306   0.9605   0.0000   0.0100  36.6918   2.3203        
+  1  1  2  67.2765  20.1739   3.3306   0.0000   0.0100   0.0000   1.1630        
+  2  1  2  74.7224  38.7524   2.1423   0.0000   0.8474   0.0000   1.3144        
+  1  2  2   0.0000   0.0000   6.0000   0.0000   0.0000   0.0000   1.0400        
+  1  2  1   0.0000   3.4110   7.7350   0.0000   0.0000   0.0000   1.0400        
+  2  2  2   0.0000  27.9213   5.8635   0.0000   0.0000   0.0000   1.0400        
+  1  1  3  61.1001  38.5283   0.4917   0.0000   4.2373   0.0000   2.1649        
+  3  1  3  72.0708  27.2877   2.3068 -10.4517   0.1000   0.0000   1.6107        
+  1  1  4  64.7353  38.2645   1.1478   0.0000   1.1834   0.0000   2.8465        
+  3  1  4  81.0672  41.9015   0.4878   0.0000   1.1019   0.0000   1.0000        
+  4  1  4  89.7621  43.0000   0.5895   0.0000   1.1155   0.0000   1.0000        
+  2  1  3  61.3200  29.1602   0.9036   0.0000   1.8063   0.0000   1.4867        
+  2  1  4  68.7361  36.7162   1.6697   0.0000   0.2000   0.0000   3.0000        
+  1  2  4   0.0000   0.0019   6.3000   0.0000   0.0000   0.0000   1.0400        
+  1  3  1  74.2533  41.5372   0.4237   0.0000   2.3660   0.0000   1.0319        
+  1  3  3  82.8809  21.0869   3.0902   0.0000   4.5310  20.1072   1.0105        
+  1  3  4  70.3730  45.0000   1.4731   0.0000   2.9000   0.0000   2.4464        
+  3  3  3  77.3022  33.4558   1.4033   0.0000   3.9048   0.0000   1.0000        
+  3  3  4  77.0669  27.6795   1.6466   0.0000   2.9000   0.0000   1.5085        
+  4  3  4  68.8583  40.4712   1.8369   0.0000   3.0072   0.0000   1.5773        
+  1  3  2  90.0000  34.3702   0.4674   0.0000   0.6101   0.0000   1.0013        
+  2  3  3  90.0000  26.3185   8.0000   0.0000   0.8615   0.0000   1.0579        
+  2  3  4  68.3253  36.1953   7.5000   0.0000   0.1000   0.0000   1.0000        
+  2  3  2  80.7909  19.0967   1.0887   0.0000   2.0594   0.0000   1.7720        
+  1  4  1  71.2077  14.1180   3.3944   0.0000   2.8702   0.0000   1.2651        
+  1  4  3  76.1064  23.7583   1.6308   0.0000   2.8701   0.0000   1.6732        
+  1  4  4  71.3624  12.8120   3.1458   0.0000   2.8701   0.0000   1.1896        
+  3  4  3  74.2922  23.0742   2.6248 -18.0069   3.0701   0.0000   1.6278        
+  3  4  4  74.0840  31.1381   1.5175  -0.9193   3.0117   0.0000   1.3541        
+  4  4  4  76.0945  32.1176   1.7767   0.0000   2.9983   0.0000   1.9677        
+  1  4  2  69.1892  14.8553   2.7174   0.0000   0.2025   0.0000   1.3071        
+  2  4  3  74.5555  45.0000   1.1948   0.0000   0.3956   0.0000   3.0000        
+  2  4  4  78.8758  45.0000   0.5964   0.0000   0.5437   0.0000   1.0000        
+  2  4  2  81.5738   7.0792   7.5000   0.0000   0.1000   0.0000   1.0000        
+  1  2  3   0.0000  21.4989   1.0000   0.0000   0.1000   0.0000   1.1358        
+  1  2  4   0.0000   0.0100   2.4974   0.0000   0.0000   0.0000   1.3777        
+  1  2  5   0.0000  15.0000   3.0000   0.0000   0.0000   0.0000   1.0400        
+  3  2  3   0.0000  15.0000   0.5640   0.0000   0.0000   0.0000   1.0400        
+  3  2  4   0.0000   1.0235   0.1000   0.0000   0.0000   0.0000   3.0000        
+  4  2  4   0.0000   0.0100   1.3170   0.0000   0.0000   0.0000   2.1165        
+  2  2  3   0.0000   2.0000   0.0839   0.0000   0.0000   0.0000   2.9374        
+  2  2  4   0.0000   0.0019   6.0000   0.0000   0.0000   0.0000   1.0400        
+  1  1  5  74.4180  33.4273   1.7018   0.1463   0.5000   0.0000   1.6178        
+  1  5  1  79.7037  28.2036   1.7073   0.1463   0.5000   0.0000   1.6453        
+  2  1  5  63.3289  29.4225   2.1326   0.0000   0.5000   0.0000   3.0000        
+  1  5  2  85.9449  38.3109   1.2492   0.0000   0.5000   0.0000   1.1000        
+  1  5  5  85.6645  40.0000   2.9274   0.1463   0.5000   0.0000   1.3830        
+  2  5  2  83.8555   5.1317   0.4377   0.0000   0.5000   0.0000   3.0000        
+  2  5  5  97.0064  32.1121   2.0242   0.0000   0.5000   0.0000   2.8568        
+  3  5  3  81.0926  30.2268   6.4132  -5.4471   2.5968   0.0000   3.0000        
+  1  5  3  70.0000  35.0000   3.4223   0.0000   1.3550   0.0000   1.2002        
+  1  3  5  57.3353  41.0012   1.0609   0.0000   1.3000   0.0000   3.0000        
+  3  3  5  83.9753  31.0715   3.5590   0.0000   0.8161   0.0000   1.1776        
+  2  3  5  89.8843  17.5000   3.3660   0.0000   2.0000   0.0000   2.0734        
+  2  6  2   0.0000  49.8261   0.2093   0.0000   2.0870   0.0000   2.2895        
+  2  2  6   0.0000  40.0366   3.1505   0.0000   1.1296   0.0000   1.1110        
+  6  2  6   0.0000   0.5047   0.8000   0.0000   0.8933   0.0000   4.6650        
+  2  6  6   0.0000   8.7037   0.0827   0.0000   3.5597   0.0000   1.1198        
+  3  6  3   0.0000   9.2317   0.1000   0.0000   1.0000   0.0000   1.0920        
+  6  3  6   0.0008  25.0000   8.0000   0.0000   1.0000   0.0000   3.0000        
+  2  3  6  66.0423   5.0000   1.0000   0.0000   1.0000   0.0000   1.2500        
+  2  6  3   0.0000   0.5000   0.1000   0.0000   1.0000   0.0000   3.0000        
+  3  3  6  70.0000  20.0000   1.0000   0.0000   1.0000   0.0000   1.2500        
+  3  7  3  90.0000  18.4167   0.6799  -8.0000   0.1310   0.0000   2.2321        
+  2  3  7  72.6004   9.6150   0.8905   0.0000   3.5473   0.0000   1.0400        
+  3  3  7  60.0000  40.0000   4.0000   0.0000   1.0000   0.0000   1.0400        
+  3  2  7   0.0000  10.0000   1.0000   0.0000   1.0000   0.0000   1.0400        
+  6  3  7  41.0995   3.2207   7.3523   0.0000   0.1101   0.0000   1.0947        
+  7  3  7  62.1312   7.5931   0.1000   0.0000   0.5154   0.0000   2.1744        
+  1  3  7  74.1394   8.5687   1.7132   0.0000  -0.6553   0.0000   2.2323        
+  2  7  3  75.0000  25.0000   2.0000   0.0000   1.0000   0.0000   1.2500        
+  3  7  7  70.0000  25.0000   2.0000   0.0000   1.0000   0.0000   1.2500        
+  3  9  3  90.0000  30.4624   2.1468   0.0000   0.0500   0.0000   1.9485        
+  9  3  9  90.0000   5.7486   5.0000   0.0000   2.0000   0.0000   1.1000        
+  3  3  9  62.9344  15.0215   4.3743   0.0000   0.6168   0.0000   1.1673        
+  3  9  9  33.7127   8.0623   3.4580   0.0000   0.0500   0.0000   2.6065        
+  2  3  9  90.0000   9.7766   8.0000   0.0000   0.0505   0.0000   1.7257        
+  1  3  9  90.0000  11.2108   1.4880   0.0000   0.5386   0.0000   2.1105        
+  3  2 10   0.0000   0.0100   0.0100   0.0000   0.0000   0.0000   1.1456        
+ 11  1 11  77.8443  49.0744   5.9913   0.0000   0.7835   0.0000   2.3020        
+  1 11  1   0.0000  19.9962   3.2299   0.0000   2.1012   0.0000   1.1537        
+  1 11 11   0.0000  25.0000   1.0000   0.0000   1.0000   0.0000   1.0400        
+ 11  1  2  69.6421  10.0000   2.0000   0.0000   1.0000   0.0000   1.0400        
+  4  9  4  90.0000  30.4624   2.1468   0.0000   0.0500   0.0000   1.9485        
+  3  9  4  90.0000  30.4624   2.1468   0.0000   0.0500   0.0000   1.9485        
+  9  4  9  90.0000   5.7486   5.0000   0.0000   2.0000   0.0000   1.1000        
+  3  4  9  62.9344  15.0215   4.3743   0.0000   0.6168   0.0000   1.1673        
+  4  3  9  62.9344  15.0215   4.3743   0.0000   0.6168   0.0000   1.1673        
+  4  4  9  62.9344  15.0215   4.3743   0.0000   0.6168   0.0000   1.1673        
+  4  4  9  62.9344  15.0215   4.3743   0.0000   0.6168   0.0000   1.1673        
+  4  9  9  33.7127   8.0623   3.4580   0.0000   0.0500   0.0000   2.6065        
+  2  4  9  90.0000   9.7766   8.0000   0.0000   0.0505   0.0000   1.7257        
+  1  4  9  90.0000  11.2108   1.4880   0.0000   0.5386   0.0000   2.1105        
+  3 12  3  90.0000   5.2360   6.0000   0.0000   1.9491   0.0000   1.0000        
+ 12  3 12  27.0723   4.9264   1.7778   0.0000   0.3851   0.0000   1.4855        
+  2  3 12   1.0000   1.5989   6.0000   0.0000   0.6668   0.0000   1.0000        
+  3  3 12  90.0000  10.0000   1.0000   0.0000   1.0000   0.0000   2.0000        
+  1  3 12  70.0000   0.0000   1.0000   0.0000   1.0000   0.0000   2.0000        
+ 73    ! Nr of torsions;at1;at2;at3;at4;;V1;V2;V3;V2(BO);vconj;n.u;n            
+  1  1  1  1   0.0734  56.5103   0.0000  -5.5375  -2.6113   0.0000   0.0000     
+  1  1  1  2   0.2335  44.4534   0.5000  -5.5234  -3.0000   0.0000   0.0000     
+  2  1  1  2  -0.0100  59.3937   0.5000  -8.0000  -3.0000   0.0000   0.0000     
+  1  1  1  3  -0.1605   9.1977  -1.0000  -2.5000  -2.1235   0.0000   0.0000     
+  2  1  1  3   0.4631  41.1330   0.7425  -7.5125  -1.1040   0.0000   0.0000     
+  3  1  1  3  -1.0000  30.3068  -1.0000  -8.5000  -0.0100   0.0000   0.0000     
+  1  1  3  1   1.0000  45.3654  -0.7715  -3.2711  -0.0100   0.0000   0.0000     
+  1  1  3  2   1.0000 150.0000   1.0000  -6.6407  -0.0100   0.0000   0.0000     
+  2  1  3  1   1.0000 143.9644   0.6834  -6.1023  -1.1889   0.0000   0.0000     
+  2  1  3  2  -1.0000  74.4134   1.0000  -3.3595  -3.5000   0.0000   0.0000     
+  1  1  3  3   1.0000 133.2711   0.5065  -2.5561  -0.0100   0.0000   0.0000     
+  2  1  3  3   0.2238  53.9090   0.8685  -2.7064  -3.5000   0.0000   0.0000     
+  3  1  3  1  -1.0000  66.2620   0.7534  -3.3533  -2.2258   0.0000   0.0000     
+  3  1  3  2   0.6489  19.7090  -1.0000  -3.3466  -0.1540   0.0000   0.0000     
+  3  1  3  3  -1.0000  86.4981   1.0000  -3.5833  -0.0610   0.0000   0.0000     
+  1  3  3  1   1.0000   0.1000   1.0000  -6.9024  -0.0100   0.0000   0.0000     
+  1  3  3  2   0.9093  23.6198  -0.1431  -2.4750  -3.5000   0.0000   0.0000     
+  2  3  3  2  -1.6122   5.0000  -1.0000  -2.5394  -0.9921   0.0000   0.0000     
+  1  3  3  3   2.5000   2.1016   0.9647  -2.6000  -0.9972   0.0000   0.0000     
+  2  3  3  3  -2.5000  75.7606  -0.6120  -7.8633  -1.2407   0.0000   0.0000     
+  3  3  3  3  -0.5000   5.0000   1.0000  -2.5000  -0.9000   0.0000   0.0000     
+  1  1  4  2   0.2700  43.1323  -0.4952  -7.6538  -1.9825   0.0000   0.0000     
+  2  1  4  2  -0.5047  82.9784   0.8701  -7.6680  -2.1051   0.0000   0.0000     
+  3  1  4  2   0.8306  15.5307   1.0000  -2.5000  -2.5261   0.0000   0.0000     
+  3  1  1  4  -0.8051  19.8307   1.0000  -3.7979  -0.9511   0.0000   0.0000     
+  4  1  1  4   1.0000  36.1913   1.0000  -3.4095  -1.7241   0.0000   0.0000     
+  1  1  4  1   1.0000  32.6616   0.3481  -6.4524  -1.6589   0.0000   0.0000     
+  3  1  4  1  -1.0000  -5.0000   1.0000  -2.5000  -1.8038   0.0000   0.0000     
+  2  1  1  4   0.7529  50.8010  -0.5000  -4.3471  -1.9000   0.0000   0.0000     
+  4  1  4  2   0.3787  13.7301   0.6579  -8.2500  -2.0202   0.0000   0.0000     
+  2  1  4  1  -1.0000  76.7186   0.1194  -8.0000  -1.5996   0.0000   0.0000     
+  0  1  2  0   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000     
+  0  2  2  0   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000     
+  0  2  3  0   0.0000   0.1000   0.0200  -2.5415   0.0000   0.0000   0.0000     
+  0  1  1  0   0.0000  60.0000   0.3000  -4.0000  -2.0000   0.0000   0.0000     
+  0  3  3  0   0.5511  25.4150   1.1330  -5.1903  -1.0000   0.0000   0.0000     
+  0  1  4  0   0.2176  40.4126   0.3535  -3.9875  -2.0051   0.0000   0.0000     
+  0  2  4  0   0.0000   0.1032   0.3000  -5.0965   0.0000   0.0000   0.0000     
+  0  3  4  0   1.1397  61.3225   0.5139  -3.8507  -2.7831   0.0000   0.0000     
+  0  4  4  0   0.7265  44.3155   1.0000  -4.4046  -2.0000   0.0000   0.0000     
+  4  1  4  4  -0.0949   8.7582   0.3310  -7.9430  -2.0000   0.0000   0.0000     
+  0  1  5  0   4.0885  78.7058   0.1174  -2.1639   0.0000   0.0000   0.0000     
+  0  5  5  0  -0.0170 -56.0786   0.6132  -2.2092   0.0000   0.0000   0.0000     
+  0  2  5  0   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000     
+  2  3  5  3   2.5000   2.5000   0.2237 -10.0000  -1.0000   0.0000   0.0000     
+  0  3  5  0  -2.5000  50.0000  -0.5000 -10.0000  -1.0000   0.0000   0.0000     
+  0  6  6  0   0.0000   0.0000   0.1200  -2.4426   0.0000   0.0000   0.0000     
+  0  2  6  0   0.0000   0.0000   0.1200  -2.4847   0.0000   0.0000   0.0000     
+  0  3  6  0   0.0000   0.0000   0.1200  -2.4703   0.0000   0.0000   0.0000     
+  1  1  1  7  -0.3232  14.3871   0.1823  -9.8682  -1.7255   0.0000   0.0000     
+  7  1  1  7  -0.1452  50.0000  -0.1915  -8.0773  -1.7255   0.0000   0.0000     
+  0  1  7  0   4.0000  45.8264   0.9000  -4.0000   0.0000   0.0000   0.0000     
+  0  7  7  0   4.0000  45.8264   0.9000  -4.0000   0.0000   0.0000   0.0000     
+  2  1  3  7  -1.5000  18.9285   0.3649  -6.1208   0.0000   0.0000   0.0000     
+  2  3  7  3   1.5000  -1.0000   0.2575  -6.2100   0.0000   0.0000   0.0000     
+  1  3  7  3  -1.4375  -0.8700   0.9861  -2.5424   0.0000   0.0000   0.0000     
+  7  3  7  3  -1.5000  21.5086  -1.0000  -4.8869   0.0000   0.0000   0.0000     
+  2  1  3  9   0.1714  69.9743   0.9170  -7.9557   0.0000   0.0000   0.0000     
+  1  1  3  9   0.2500  76.5218   1.0000  -2.5503   0.0000   0.0000   0.0000     
+  3  1  3  9  -0.2500  50.5929  -0.2500  -6.9285   0.0000   0.0000   0.0000     
+  2  3  9  3  -0.2500   0.0100  -0.5000  -4.6984   0.0000   0.0000   0.0000     
+  1  1  1 11   0.5000   0.1000   0.4683 -11.5274  -1.7255   0.0000   0.0000     
+  2  1  1 11   0.0000  49.3871   0.2000 -10.5765  -1.7255   0.0000   0.0000     
+ 11  1  1 11  -0.5000  95.4727  -0.2080  -4.8579  -1.7255   0.0000   0.0000     
+  0  1 11  0   4.0000  45.8264   0.9000  -4.0000   0.0000   0.0000   0.0000     
+  0 11 11  0   4.0000  45.8264   0.8897  -4.0000   0.0000   0.0000   0.0000     
+  2  1  4  9   0.1714  69.9743   0.9170  -7.9557   0.0000   0.0000   0.0000     
+  1  1  4  9   0.2500  76.5218   1.0000  -2.5503   0.0000   0.0000   0.0000     
+  3  1  4  9  -0.2500  50.5929  -0.2500  -6.9285   0.0000   0.0000   0.0000     
+  4  1  4  9  -0.2500  50.5929  -0.2500  -6.9285   0.0000   0.0000   0.0000     
+  4  1  3  9  -0.2500  50.5929  -0.2500  -6.9285   0.0000   0.0000   0.0000     
+  2  4  9  3  -0.2500   0.0100  -0.5000  -4.6984   0.0000   0.0000   0.0000     
+  2  4  9  4  -0.2500   0.0100  -0.5000  -4.6984   0.0000   0.0000   0.0000     
+  0    ! Nr of hydrogen bonds;at1;at2;at3;Rhb;Dehb;vhb1                         
diff --git a/data/benchmarks/water/ffield_acks2.water b/data/benchmarks/water/ffield_acks2.water
new file mode 100644
index 0000000000000000000000000000000000000000..ddf2d3cbf5f8c76d3203ee4a326ee8db23fb246e
--- /dev/null
+++ b/data/benchmarks/water/ffield_acks2.water
@@ -0,0 +1,364 @@
+Reactive MD-force field: Water                                                  
+ 39       ! Number of general parameters                                        
+   50.0000 !Overcoordination parameter                                          
+    9.5469 !Overcoordination parameter                                          
+   26.5405 !Valency angle conjugation parameter                                 
+    1.7224 !Triple bond stabilisation parameter                                 
+    6.8702 !Triple bond stabilisation parameter                                 
+   60.4850 !C2-correction                                                       
+    1.0588 !Undercoordination parameter                                         
+    4.6000 !Triple bond stabilisation parameter                                 
+   12.1176 !Undercoordination parameter                                         
+   13.3056 !Undercoordination parameter                                         
+  -70.5044 !Triple bond stabilization energy                                    
+    0.0000 !Lower Taper-radius                                                  
+   10.0000 !Upper Taper-radius                                                  
+    2.8793 !Not used                                                            
+   33.8667 !Valency undercoordination                                           
+    6.0891 !Valency angle/lone pair parameter                                   
+    1.0563 !Valency angle                                                       
+    2.0384 !Valency angle parameter                                             
+    6.1431 !Not used                                                            
+    6.9290 !Double bond/angle parameter                                         
+    0.3989 !Double bond/angle parameter: overcoord                              
+    3.9954 !Double bond/angle parameter: overcoord                              
+   -2.4837 !Not used                                                            
+    5.7796 !Torsion/BO parameter                                                
+   10.0000 !Torsion overcoordination                                            
+    1.9487 !Torsion overcoordination                                            
+   -1.2327 !Conjugation 0 (not used)                                            
+    2.1645 !Conjugation                                                         
+    1.5591 !vdWaals shielding                                                   
+    0.1000 !Cutoff for bond order (*100)                                        
+    2.1365 !Valency angle conjugation parameter                                 
+    0.6991 !Overcoordination parameter                                          
+   50.0000 !Overcoordination parameter                                          
+    1.8512 !Valency/lone pair parameter                                         
+  548.6451 !Softness                                                            
+   20.0000 !Not used                                                            
+    5.0000 !Molecular energy (not used)                                         
+    0.0000 !Molecular energy (not used)                                         
+    2.6962 !Valency angle conjugation parameter                                 
+ 15    ! Nr of atoms; cov.r; valency;a.m;Rvdw;Evdw;gammaEEM;cov.r2;#            
+            alfa;gammavdW;valency;Eunder;Eover;chiEEM;etaEEM;n.u.               
+            cov r3;Elp;Heat inc.;n.u.;n.u.;n.u.;n.u.                            
+            ov/un;val1;n.u.;val3,vval4                                          
+ C    1.3817   4.0000  12.0000   1.8903   0.1838   0.9000   1.1341   4.0000     
+      9.7559   2.1346   4.0000  34.9350  79.5548   5.9666   7.0000   0.0000     
+      1.2114   0.0000 202.5551   8.9539  34.9289  13.5366   0.8563   0.0000     
+     -2.8983   2.5000   1.0564   4.0000   2.9663   0.0000   0.0000   0.0000     
+ H    0.8930   1.0000   1.0080   1.3550   0.0930   0.8203  -0.1000   1.0000     
+      8.2230  33.2894   1.0000   0.0000 121.1250   3.7248   9.6093   1.0000     
+     -0.1000   0.0000  61.6606   3.0408   2.4197   0.0003   3.4114   0.0000     
+    -19.4571   4.2733   1.0338   1.0000   2.8793   0.0000   0.0000   0.0000     
+ O    1.2450   2.0000  15.9990   2.3890   0.1000   1.0898   1.0548   6.0000     
+      9.7300  13.8449   4.0000  37.5000 116.0768   8.5000   8.3122   2.0000     
+      0.9049   0.4056  59.0626   3.5027   0.7640   0.0021   0.9745   0.0000     
+     -3.5500   2.9000   1.0493   4.0000   2.9225   0.0000   0.0000   0.0000     
+ N    1.2333   3.0000  14.0000   1.9324   0.1376   0.8596   1.1748   5.0000     
+     10.0667   7.8431   4.0000  32.2482 100.0000   6.8418   6.3404   2.0000     
+      1.0433  13.7673 119.9837   2.1961   3.0696   2.7683   0.9745   0.0000     
+     -4.3875   2.6192   1.0183   4.0000   2.8793   0.0000   0.0000   0.0000     
+ S    1.9405   2.0000  32.0600   2.0677   0.2099   1.0336   1.5479   6.0000     
+      9.9575   4.9055   4.0000  52.9998 112.1416   6.5000   8.2545   2.0000     
+      1.4601   9.7177  71.1843   5.7487  23.2859  12.7147   0.9745   0.0000     
+    -11.0000   2.7466   1.0338   6.2998   2.8793   0.0000   0.0000   0.0000     
+ Si   2.0276   4.0000  28.0600   2.2042   0.1322   0.8218   1.5758   4.0000     
+     11.9413   2.0618   4.0000  11.8211 136.4845   1.8038   7.3852   0.0000     
+     -1.0000   0.0000 126.5182   6.4918   8.5961   0.2368   0.8563   0.0000     
+     -3.8112   3.1873   1.0338   6.2998   2.5791   0.0000   0.0000   0.0000     
+ Pt   1.9907   3.0000 195.0800   1.9980   0.2452   0.8218  -1.0000   3.0000     
+     12.8669   3.2118   3.0000   0.0000   0.0000   1.8038   7.3852   0.0000     
+     -1.0000   0.0000 142.6300   6.2293   5.2294   0.1542   0.8563   0.0000     
+     -6.7740   2.9867   1.0338   6.2998   2.5791   0.0000   0.0000   0.0000     
+ Zr   2.1000   4.0000  91.2240   2.1970   0.2542   0.8218  -1.0000   4.0000     
+     12.8545   3.5938   4.0000   0.0000   0.0000   1.8038   7.3852   0.0000     
+     -1.0000   0.0000 107.6300   6.2293   5.2294   0.1542   0.8563   0.0000     
+     -3.2224   2.9867   1.0338   6.2998   2.5791   0.0000   0.0000   0.0000     
+ Ni   1.8503   2.0000  58.6900   1.9219   0.1582   0.8218  -1.0000   2.0000     
+     12.1238   4.0351   2.0000   0.0000   0.0000   1.8038   7.3852   0.0000     
+     -1.0000   0.0000  95.6300   6.2293   5.2294   0.1542   0.8563   0.0000     
+     -3.2224   2.9867   1.0338   6.2998   2.5791   0.0000   0.0000   0.0000     
+ Au   1.8503   1.0000 196.9665   1.9219   0.1582   0.8218  -1.0000   1.0000     
+     12.1238   4.0351   1.0000   0.0000   0.0000   1.8038   7.3852   0.0000     
+     -1.0000   0.0000  72.6300   6.2293   5.2294   0.1542   0.8563   0.0000     
+     -3.2224   2.9867   1.0338   6.2998   2.5791   0.0000   0.0000   0.0000     
+ V    2.2657   3.0000  50.9415   1.7992   0.3005   0.6743   0.1000   5.0000     
+     12.3879   5.2243   3.0000   0.0000   0.0000  -0.3628   6.6023   0.0000     
+     -1.0000   0.0000 117.6300  23.1946   6.5795   0.0000   0.8563   0.0000     
+     -3.5389   1.5012   1.0338   3.0000   3.6411   0.0000   0.0000   0.0000     
+ Bi   2.1949   3.0000 208.9804   2.4429   0.1607   0.4960   0.0535   5.0000     
+     12.9571  35.5167   3.0000   0.0000   0.0000  -0.1926   6.4153   0.0000     
+     -1.0000   0.5785  52.6300   3.8978   0.9856   0.0314   0.8563   0.0000     
+     -2.5000   5.0597   1.0338   6.0000   2.5791   0.0000   0.0000   0.0000     
+ Ti   0.1000   4.0000  47.8800   2.0000   0.1659   0.6037   0.1000   4.0000     
+     13.2535   4.0063   4.0000  -5.0000   0.0000  -0.1864   5.9304   0.0000     
+     -1.0000   0.0000 129.6300  22.8461   1.8515   0.0064   0.8563   0.0000     
+     -3.4122   3.2711   1.0338   6.2998   2.2632   0.0000   0.0000   0.0000     
+ Mo   2.4710   5.6504  95.9400   1.8000   0.3285   1.0000   0.1000   6.0000     
+     13.0000  45.0000   4.0000   0.0000   0.0000   0.6062   6.1484   0.0000     
+      0.1000   0.0000 152.6300   3.7659   0.0689   2.9902   0.8563   0.0000     
+    -16.7660   3.1072   1.0338   8.0000   3.4590   0.0000   0.0000   0.0000     
+ X   -0.1000   2.0000   1.0080   2.0000   0.0000   1.0000  -0.1000   6.0000     
+     10.0000   2.5000   4.0000   0.0000   0.0000   8.5000   1.5000   0.0000     
+     -0.1000   0.0000  -2.3700   8.7410  13.3640   0.6690   0.9745   0.0000     
+    -11.0000   2.7466   1.0338   2.0000   2.8793   0.0000   0.0000   0.0000     
+ 40      ! Nr of bonds; Edis1;LPpen;n.u.;pbe1;pbo5;13corr;pbo6                  
+                         pbe2;pbo3;pbo4;n.u.;pbo1;pbo2;ovcorr                   
+  1  1 158.2004  99.1897  78.0000  -0.7738  -0.4550   1.0000  37.6117   0.4147  
+         0.4590  -0.1000   9.1628   1.0000  -0.0777   6.7268   1.0000   0.0000  
+  1  2 169.4760   0.0000   0.0000  -0.6083   0.0000   1.0000   6.0000   0.7652  
+         5.2290   1.0000   0.0000   1.0000  -0.0500   6.9136   0.0000   0.0000  
+  2  2 153.3934   0.0000   0.0000  -0.4600   0.0000   1.0000   6.0000   0.7300  
+         6.2500   1.0000   0.0000   1.0000  -0.0790   6.0552   0.0000   0.0000  
+  1  3 158.6946 107.4583  23.3136  -0.4240  -0.1743   1.0000  10.8209   1.0000  
+         0.5322  -0.3113   7.0000   1.0000  -0.1447   5.2450   0.0000   0.0000  
+  3  3 142.2858 145.0000  50.8293   0.2506  -0.1000   1.0000  29.7503   0.6051  
+         0.3451  -0.1055   9.0000   1.0000  -0.1225   5.5000   1.0000   0.0000  
+  1  4 134.1215 140.2179  79.9745   0.0163  -0.1428   1.0000  27.0617   0.2000  
+         0.1387  -0.3681   7.1611   1.0000  -0.1000   5.0825   1.0000   0.0000  
+  3  4 130.8596 169.4551  40.0000   0.3837  -0.1639   1.0000  35.0000   0.2000  
+         1.0000  -0.3579   7.0004   1.0000  -0.1193   6.8773   1.0000   0.0000  
+  4  4 157.9384  82.5526 152.5336   0.4010  -0.1034   1.0000  12.4261   0.5828  
+         0.1578  -0.1509  11.9186   1.0000  -0.0861   5.4271   1.0000   0.0000  
+  2  3 160.0000   0.0000   0.0000  -0.5725   0.0000   1.0000   6.0000   0.5626  
+         1.1150   1.0000   0.0000   0.0000  -0.0920   4.2790   0.0000   0.0000  
+  2  4 231.8173   0.0000   0.0000  -0.3364   0.0000   1.0000   6.0000   0.4402  
+         8.8910   1.0000   0.0000   1.0000  -0.0327   6.5754   0.0000   0.0000  
+  1  5 128.9942  74.5848  55.2528   0.1035  -0.5211   1.0000  18.9617   0.6000  
+         0.2949  -0.2398   8.1175   1.0000  -0.1029   5.6731   1.0000   0.0000  
+  2  5 151.5159   0.0000   0.0000  -0.4721   0.0000   1.0000   6.0000   0.6000  
+         9.4366   1.0000   0.0000   1.0000  -0.0290   7.0050   1.0000   0.0000  
+  3  5   0.0000   0.0000   0.0000   0.5563  -0.4038   1.0000  49.5611   0.6000  
+         0.4259  -0.4577  12.7569   1.0000  -0.1100   7.1145   1.0000   0.0000  
+  4  5   0.0000   0.0000   0.0000   0.4438  -0.2034   1.0000  40.3399   0.6000  
+         0.3296  -0.3153   9.1227   1.0000  -0.1805   5.6864   1.0000   0.0000  
+  5  5  96.1871  93.7006  68.6860   0.0955  -0.4781   1.0000  17.8574   0.6000  
+         0.2723  -0.2373   9.7875   1.0000  -0.0950   6.4757   1.0000   0.0000  
+  6  6 109.1904  70.8314  30.0000   0.2765  -0.3000   1.0000  16.0000   0.1583  
+         0.2804  -0.1994   8.1117   1.0000  -0.0675   8.2993   0.0000   0.0000  
+  2  6 137.1002   0.0000   0.0000  -0.1902   0.0000   1.0000   6.0000   0.4256  
+        17.7186   1.0000   0.0000   1.0000  -0.0377   6.4281   0.0000   0.0000  
+  3  6 191.1743  52.0733  43.3991  -0.2584  -0.3000   1.0000  36.0000   0.8764  
+         1.0248  -0.3658   4.2151   1.0000  -0.5004   4.2605   1.0000   0.0000  
+  4  6 185.4488  39.2832  43.3991  -0.1922  -0.3000   1.0000  36.0000   0.8217  
+         0.8538  -0.3887   4.4334   1.0000  -0.5241   4.4529   1.0000   0.0000  
+  7  7  90.1462   0.0000   0.0000   0.0004  -0.2000   0.0000  16.0000   0.3484  
+         1.0000  -0.2000  15.0000   1.0000  -0.1014   5.7631   0.0000   0.0000  
+  8  8  85.2900   0.0000   0.0000   0.0004  -0.2000   0.0000  16.0000   0.5438  
+         1.0000  -0.2000  15.0000   1.0000  -0.1001   5.5699   0.0000   0.0000  
+  9  9  73.6182   0.0000   0.0000   0.0004  -0.2000   0.0000  16.0000   0.3418  
+         1.0000  -0.2000  15.0000   1.0000  -0.1015   5.7850   0.0000   0.0000  
+ 10 10  73.6182   0.0000   0.0000   0.0004  -0.2000   0.0000  16.0000   0.3418  
+         1.0000  -0.2000  15.0000   1.0000  -0.1015   5.7850   0.0000   0.0000  
+ 11 11  36.2751   0.0000   0.0000   0.8059  -0.3000   0.0000  16.0000   0.1826  
+         0.3414  -0.3000  16.0000   1.0000  -0.0717   7.9108   0.0000   0.0000  
+  3 11 106.8008  67.5543   0.0000   0.0323  -0.3000   1.0000  36.0000   0.1000  
+         0.2670  -0.3402  16.0000   1.0000  -0.1761   4.6698   1.0000   0.0000  
+  2 11   0.0000   0.0000   0.0000  -0.2872  -0.3000   1.0000  36.0000   0.0082  
+         1.7973  -0.2500  20.0000   1.0000  -0.2578   6.5219   1.0000   0.0000  
+  1 11   0.0000   0.0000   0.0000  -0.2872  -0.3000   1.0000  36.0000   0.0082  
+         1.7973  -0.2500  20.0000   1.0000  -0.2578   6.5219   1.0000   0.0000  
+ 12 12  66.0677   0.0000   0.0000  -0.9557  -0.2000   0.0000  16.0000   0.2865  
+         0.5847  -0.2000  15.0000   1.0000  -0.0856   5.2857   0.0000   0.0000  
+  3 12 152.2407  57.6204   0.0000  -0.8033  -0.3000   1.0000  36.0000   0.0498  
+         1.8097  -0.3800  16.0000   1.0000  -0.2379   8.0000   1.0000   0.0000  
+  2 12  95.9209   0.0000   0.0000  -0.0153  -0.3000   1.0000  36.0000   0.0100  
+         1.0000  -0.2062   8.6647   1.0000  -0.1911   4.0000   1.0000   0.0000  
+  1 12  78.9091  40.6322   0.0000   0.0040  -0.3000   1.0000  36.0000   0.0384  
+         0.0904  -0.1209  12.3682   1.0000  -0.1613   4.3849   1.0000   0.0000  
+ 13 13  71.3016  10.0000   0.0000  -0.1571  -0.2000   0.0000  16.0000   0.3311  
+         0.1822  -0.2000  15.0000   1.0000  -0.1860   6.5172   0.0000   0.0000  
+  3 13 112.7130  29.8084   0.0000  -0.9010  -0.3000   1.0000  36.0000   0.5508  
+         0.1006  -0.2492  16.9476   1.0000  -0.1919   5.4797   1.0000   0.0000  
+  1 13   0.0000   0.0000   0.0000  -0.2872  -0.3000   1.0000  36.0000   0.0082  
+         1.7973  -0.2500  20.0000   1.0000  -0.2578   6.5219   1.0000   0.0000  
+  2 13   0.0000   0.0000   0.0000  -0.2872  -0.3000   1.0000  36.0000   0.0082  
+         1.7973  -0.2500  20.0000   1.0000  -0.2578   6.5219   1.0000   0.0000  
+  1 14   0.5356   0.9614   0.0000   0.3817  -0.3000   1.0000  36.0000   0.2142  
+         0.6116  -0.2579   6.1366   1.0000  -0.0913   6.6008   1.0000   0.0000  
+  2 14   0.0000   0.0000   0.0000  -0.2872  -0.3000   1.0000  36.0000   0.0082  
+         1.7973  -0.3027   4.6243   1.0000  -0.4578   3.5219   1.0000   0.0000  
+  3 14 112.7070  10.0000 135.5011   0.9277  -0.2354   1.0000  19.1731   1.2334  
+         0.9822  -0.1837   7.2216   1.0000  -0.1264   6.1257   1.0000   0.0000  
+ 14 14  44.6382   0.0000   0.0000   1.0000  -0.3000   0.0000  16.0000   0.2890  
+         0.3384  -0.3000  16.0000   1.0000  -0.1862   7.4588   0.0000   0.0000  
+ 12 14  50.0000   0.0000   0.0000   0.1000  -0.3000   0.0000  16.0000   0.3000  
+         1.0000  -0.3000  16.0000   1.0000  -0.2000   8.0000   0.0000   0.0000  
+ 20    ! Nr of off-diagonal terms; Ediss;Ro;gamma;rsigma;rpi;rpi2               
+  1  2   0.1239   1.4004   9.8467   1.1210  -1.0000  -1.0000                    
+  2  3   0.0283   1.2885  10.9190   0.9215  -1.0000  -1.0000                    
+  2  4   0.1059   1.8290   9.7818   0.9598  -1.0000  -1.0000                    
+  1  3   0.1156   1.8520   9.8317   1.2854   1.1352   1.0706                    
+  1  4   0.1447   1.8766   9.7990   1.3436   1.1885   1.1363                    
+  3  4   0.1048   2.0003  10.1220   1.3173   1.1096   1.0206                    
+  2  6   0.0470   1.6738  11.6877   1.1931  -1.0000  -1.0000                    
+  3  6   0.1263   1.8163  10.6833   1.6266   1.2052  -1.0000                    
+  1 11   0.1995   2.2133  13.0000   0.0102   1.4868  -1.0000                    
+  2 11   0.1319   1.5855  12.5457   0.0099   1.5065  -1.0000                    
+  3 11   0.0813   1.8649  10.8791   1.6498   1.6445  -1.0000                    
+  1 12   0.4235   1.7716  11.3664   1.8000   1.7212  -1.0000                    
+  2 12   0.0754   1.6033  12.4204   1.6896  -1.5000  -1.0000                    
+  3 12   0.1648   2.1260  11.2425   2.0692   1.6939  -1.0000                    
+  2 13   0.1340   1.8546  11.5784   1.0000  -1.0000  -1.0000                    
+  3 13   0.1280   1.8000  10.5743   1.7358   1.5296  -1.0000                    
+  1 13   0.1301   1.9382  11.1255   0.0100  -1.0000  -1.0000                    
+  1 14   0.1495   2.0794  12.2376   0.0100   1.4060  -1.0000                    
+  2 14   0.0795   1.6794  11.2376   0.0100   1.2060  -1.0000                    
+  3 14   0.2101   2.0342  10.4729   1.6019   1.4781   1.6548                    
+ 97    ! Nr of angles;at1;at2;at3;Thetao,o;ka;kb;pv1;pv2                        
+  1  1  1  59.0573  30.7029   0.7606   0.0000   0.7180   6.2933   1.1244        
+  1  1  2  65.7758  14.5234   6.2481   0.0000   0.5665   0.0000   1.6255        
+  2  1  2  70.2607  25.2202   3.7312   0.0000   0.0050   0.0000   2.7500        
+  1  2  2   0.0000   0.0000   6.0000   0.0000   0.0000   0.0000   1.0400        
+  1  2  1   0.0000   3.4110   7.7350   0.0000   0.0000   0.0000   1.0400        
+  2  2  2   0.0000  27.9213   5.8635   0.0000   0.0000   0.0000   1.0400        
+  1  1  3  49.6811   7.1713   4.3889   0.0000   0.7171  10.2661   1.0463        
+  3  1  3  77.7473  40.1718   2.9802 -25.3063   1.6170 -46.1315   2.2503        
+  1  1  4  66.1305  12.4661   7.0000   0.0000   3.0000  50.0000   1.1880        
+  3  1  4  73.9544  12.4661   7.0000   0.0000   3.0000   0.0000   1.1880        
+  4  1  4  64.1581  12.4661   7.0000   0.0000   3.0000   0.0000   1.1880        
+  2  1  3  65.0000  13.8815   5.0583   0.0000   0.4985   0.0000   1.4900        
+  2  1  4  74.2929  31.0883   2.6184   0.0000   0.0755   0.0000   1.0500        
+  1  2  4   0.0000   0.0019   6.3000   0.0000   0.0000   0.0000   1.0400        
+  1  3  1  73.5312  44.7275   0.7354   0.0000   3.0000   0.0000   1.0684        
+  1  3  3  79.4761  36.3701   1.8943   0.0000   0.7351  67.6777   3.0000        
+  1  3  4  82.4890  31.4554   0.9953   0.0000   1.6310   0.0000   1.0783        
+  3  3  3  80.7324  30.4554   0.9953   0.0000   1.6310  50.0000   1.0783        
+  3  3  4  84.3637  31.4554   0.9953   0.0000   1.6310   0.0000   1.0783        
+  4  3  4  89.7071  31.4554   0.9953   0.0000   1.6310   0.0000   1.1519        
+  1  3  2  70.1880  20.9562   0.3864   0.0000   0.0050   0.0000   1.6924        
+  2  3  3  75.6935  50.0000   2.0000   0.0000   1.0000   0.0000   1.1680        
+  2  3  4  75.6201  18.7919   0.9833   0.0000   0.1218   0.0000   1.0500        
+  2  3  2  85.8000   9.8453   2.2720   0.0000   2.8635   0.0000   1.5800        
+  1  4  1  66.0330  22.0295   1.4442   0.0000   1.6777   0.0000   1.0500        
+  1  4  3 103.3204  33.0381   0.5787   0.0000   1.6777   0.0000   1.0500        
+  1  4  4 104.1335   8.6043   1.6495   0.0000   1.6777   0.0000   1.0500        
+  3  4  3  74.1978  42.1786   1.7845 -18.0069   1.6777   0.0000   1.0500        
+  3  4  4  74.8600  43.7354   1.1572  -0.9193   1.6777   0.0000   1.0500        
+  4  4  4  75.0538  14.8267   5.2794   0.0000   1.6777   0.0000   1.0500        
+  1  4  2  69.1106  25.5067   1.1003   0.0000   0.0222   0.0000   1.0369        
+  2  4  3  81.3686  40.0712   2.2396   0.0000   0.0222   0.0000   1.0369        
+  2  4  4  83.0104  43.4766   1.5328   0.0000   0.0222   0.0000   1.0500        
+  2  4  2  70.8687  12.0168   5.0132   0.0000   0.0222   0.0000   1.1243        
+  1  2  3   0.0000  25.0000   3.0000   0.0000   1.0000   0.0000   1.0400        
+  1  2  4   0.0000   0.0019   6.0000   0.0000   0.0000   0.0000   1.0400        
+  1  2  5   0.0000   0.0019   6.0000   0.0000   0.0000   0.0000   1.0400        
+  3  2  3   0.0000  15.0000   2.8900   0.0000   0.0000   0.0000   2.8774        
+  3  2  4   0.0000   0.0019   6.0000   0.0000   0.0000   0.0000   1.0400        
+  4  2  4   0.0000   0.0019   6.0000   0.0000   0.0000   0.0000   1.0400        
+  2  2  3   0.0000   8.5744   3.0000   0.0000   0.0000   0.0000   1.0421        
+  2  2  4   0.0000   0.0019   6.0000   0.0000   0.0000   0.0000   1.0400        
+  1  1  5  74.9397  25.0560   1.8787   0.1463   0.0559   0.0000   1.0400        
+  1  5  1  86.9521  36.9951   2.0903   0.1463   0.0559   0.0000   1.0400        
+  2  1  5  74.9397  25.0560   1.8787   0.0000   0.0000   0.0000   1.0400        
+  1  5  2  86.1791  36.9951   2.0903   0.0000   0.0000   0.0000   1.0400        
+  1  5  5  85.3644  36.9951   2.0903   0.1463   0.0559   0.0000   1.0400        
+  2  5  2  93.1959  36.9951   2.0903   0.0000   0.0000   0.0000   1.0400        
+  2  5  5  84.3331  36.9951   2.0903   0.0000   0.0000   0.0000   1.0400        
+  6  6  6  69.3456  21.7361   1.4283   0.0000  -0.2101   0.0000   1.3241        
+  2  6  6  75.6168  21.5317   1.0435   0.0000   2.5179   0.0000   1.0400        
+  2  6  2  78.3939  20.9772   0.8630   0.0000   2.8421   0.0000   1.0400        
+  3  6  6  70.3016  15.4081   1.3267   0.0000   2.1459   0.0000   1.0400        
+  2  6  3  73.8232  16.6592   3.7425   0.0000   0.8613   0.0000   1.0400        
+  3  6  3  90.0344   7.7656   1.7264   0.0000   0.7689   0.0000   1.0400        
+  6  3  6  22.1715   3.6615   0.3160   0.0000   4.1125   0.0000   1.0400        
+  2  3  6  83.7634   5.6693   2.7780   0.0000   1.6982   0.0000   1.0400        
+  3  3  6  73.4663  25.0761   0.9143   0.0000   2.2466   0.0000   1.0400        
+  2  2  6   0.0000  47.1300   6.0000   0.0000   1.6371   0.0000   1.0400        
+  6  2  6   0.0000  31.5209   6.0000   0.0000   1.6371   0.0000   1.0400        
+  3  2  6   0.0000  31.0427   4.5625   0.0000   1.6371   0.0000   1.0400        
+  2  2  5   0.0000   0.0019   6.0000   0.0000   0.0000   0.0000   1.0400        
+  3 11  3  62.4906  31.5023   1.3328   0.0000   2.8731   0.0000   1.0794        
+ 11  3 11  31.0790  19.3435   0.4919   0.0000   2.9625   0.0000   3.0000        
+  3  3 11 100.0000  14.7642   7.0000   0.0000   1.0585   0.0000   1.1599        
+  1  3 11  60.7895  13.6681   0.7546   0.0000   2.1747   0.0000   2.9508        
+  2  3 11 100.0000   5.0000   1.4335   0.0000   1.2363   0.0000   5.0000        
+  3 12 12  23.8296   8.9089   7.0000   0.0000   1.0000   0.0000   2.8891        
+  3 12  3  87.0764  19.4489   2.5080   0.0000   2.6056   0.0000   3.0000        
+ 12  3 12  72.7369  13.7522   5.0243   0.0000   2.9700   0.0000   1.5506        
+  3  3 12  68.8771  10.5000   2.5500   0.0000   2.5729   0.0000   1.5892        
+  2  3 12  99.5836   5.4142   2.2105   0.0000   1.0513   0.0000   1.1000        
+  1  3 12  90.0000  12.1772   2.2055   0.0000   1.9064   0.0000   2.6056        
+  1  1 12  71.1708  32.6379   0.4516   0.0000   2.1609   0.0000   1.1000        
+  1 12  3  90.0000  45.0000   0.9335   0.0000   0.2140   0.0000   1.4846        
+  1 12  1  87.6204  45.0000   1.2740   0.0000   1.1519   0.0000   1.1000        
+  3  1 12  54.7020   3.2967   7.0000   0.0000   2.0408   0.0000   2.4032        
+  2 12  3  90.0000  28.2099   1.8036   0.0000   1.5461   0.0000   1.2304        
+  2 12  2  90.0000  36.3001   0.6409   0.0000   3.0000   0.0000   1.7755        
+  1 12  2  89.5835  45.0000   0.8465   0.0000   1.2118   0.0000   2.2282        
+  2  1 12  68.7714  22.9669   0.4631   0.0000   2.4269   0.0000   1.4680        
+  2  2 12   0.0000  30.2898   3.9181   0.0000   0.9914   0.0000   1.3121        
+  3  2 12   0.0000   1.0000   4.1706   0.0000   1.0100   0.0000   1.1000        
+  1  2 12   0.0000   1.0000   3.9722   0.0000   1.0075   0.0000   1.2984        
+  3 13  3  73.6321  10.6453   2.7693   0.0000   0.0500   0.0000   1.9906        
+ 13  3 13 100.0000   5.0270   5.0000   0.0000   1.2768   0.0000   2.0630        
+  3  3 13  52.3127  40.0000   1.1362   0.0000   1.5100   0.0000   1.1000        
+  3 13 13  66.6695   0.0036   3.2646   0.0000   0.0581   0.0000   1.3741        
+  2  3 13 100.0000   3.8927   8.0000   0.0000   2.0000   0.0000   1.1000        
+  1  3 13  96.6040   9.4537   8.0000   0.0000   0.3285   0.0000   4.0000        
+  3 14  3  79.6765  50.0000   1.0502  -0.0016   0.1000   0.0000   1.4583        
+ 14  3 14  20.2100  37.6165   0.6059   0.0000   0.1531   0.0000   2.0586        
+  3  3 14  38.5570  11.9307   0.9911   0.0000   0.8422   0.0000   1.0500        
+  3 14 14   5.8342   0.0724   0.1000   0.0000   0.5490   0.0000   1.7839        
+  2  3 14  81.8943   7.2820   2.1490   0.0000   0.6873   0.0000   3.2184        
+  1  3 14  75.5634   8.3289   1.0236   0.0000   2.0875   0.0000   1.0500        
+ 12  3 14  30.0000   5.0000   0.5000   0.0000   0.5000   0.0000   1.2500        
+ 47    ! Nr of torsions;at1;at2;at3;at4;;V1;V2;V3;V2(BO);vconj;n.u;n            
+  1  1  1  1  -0.2500  34.7453   0.0288  -6.3507  -1.6000   0.0000   0.0000     
+  1  1  1  2  -0.2500  29.2131   0.2945  -4.9581  -2.1802   0.0000   0.0000     
+  2  1  1  2  -0.2500  31.2081   0.4539  -4.8923  -2.2677   0.0000   0.0000     
+  1  1  1  3  -0.3495  22.2142  -0.2959  -2.5000  -1.9066   0.0000   0.0000     
+  2  1  1  3   0.0646  24.3195   0.6259  -3.9603  -1.0000   0.0000   0.0000     
+  3  1  1  3  -0.5456   5.5756   0.8433  -5.1924  -1.0180   0.0000   0.0000     
+  1  1  3  1   1.7555  27.9267   0.0072  -2.6533  -1.0000   0.0000   0.0000     
+  1  1  3  2  -1.4358  36.7830  -1.0000  -8.1821  -1.0000   0.0000   0.0000     
+  2  1  3  1  -1.3959  34.5053   0.7200  -2.5714  -2.1641   0.0000   0.0000     
+  2  1  3  2  -2.5000  70.0597   1.0000  -3.5539  -2.9929   0.0000   0.0000     
+  1  1  3  3   0.6852  11.2819  -0.4784  -2.5000  -2.1085   0.0000   0.0000     
+  2  1  3  3   0.1933  80.0000   1.0000  -4.0590  -3.0000   0.0000   0.0000     
+  3  1  3  1  -1.9889  76.4820  -0.1796  -3.8301  -3.0000   0.0000   0.0000     
+  3  1  3  2   0.2160  72.7707  -0.7087  -4.2100  -3.0000   0.0000   0.0000     
+  3  1  3  3  -2.5000  71.0772   0.2542  -3.1631  -3.0000   0.0000   0.0000     
+  1  3  3  1   2.5000  -0.6002   1.0000  -3.4297  -2.8858   0.0000   0.0000     
+  1  3  3  2  -2.5000  -3.3822   0.7004  -5.4467  -2.9586   0.0000   0.0000     
+  2  3  3  2   2.5000  -4.0000   0.9000  -2.5000  -1.0000   0.0000   0.0000     
+  1  3  3  3   1.2329  -4.0000   1.0000  -2.5000  -1.7479   0.0000   0.0000     
+  2  3  3  3   0.8302  -4.0000  -0.7763  -2.5000  -1.0000   0.0000   0.0000     
+  3  3  3  3  -2.5000  -4.0000   1.0000  -2.5000  -1.0000   0.0000   0.0000     
+  0  1  2  0   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000     
+  0  2  2  0   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000     
+  0  2  3  0   0.0000   0.1000   0.0200  -2.5415   0.0000   0.0000   0.0000     
+  0  1  1  0   0.0000  50.0000   0.3000  -4.0000  -2.0000   0.0000   0.0000     
+  0  3  3  0   0.5511  25.4150   1.1330  -5.1903  -1.0000   0.0000   0.0000     
+  0  1  4  0  -2.4242 128.1636   0.3739  -6.6098  -2.0000   0.0000   0.0000     
+  0  2  4  0   0.0000   0.1000   0.0200  -2.5415   0.0000   0.0000   0.0000     
+  0  3  4  0   1.4816  55.6641   0.0004  -7.0465  -2.7203   0.0000   0.0000     
+  0  4  4  0  -0.3244  27.7086   0.0039  -2.8272  -2.0000   0.0000   0.0000     
+  4  1  4  4  -5.5181   8.9706   0.0004  -6.1782  -2.0000   0.0000   0.0000     
+  0  1  5  0   3.3423  30.3435   0.0365  -2.7171   0.0000   0.0000   0.0000     
+  0  5  5  0  -0.0555 -42.7738   0.1515  -2.2056   0.0000   0.0000   0.0000     
+  0  2  5  0   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000     
+  0  6  6  0   0.0000   0.0000   0.1200  -2.4426   0.0000   0.0000   0.0000     
+  0  2  6  0   0.0000   0.0000   0.1200  -2.4847   0.0000   0.0000   0.0000     
+  0  3  6  0   0.0000   0.0000   0.1200  -2.4703   0.0000   0.0000   0.0000     
+  2  1  3 14   1.6297  56.8132   0.3398  -2.6912  -2.1000   0.0000   0.0000     
+  1  1  3 14  -0.0427  13.4096   0.9351  -6.5245  -2.1000   0.0000   0.0000     
+  2  3 14  3   2.5000  11.6208   1.0000  -9.0000  -1.0000   0.0000   0.0000     
+  2  1  3 12  -0.2500  45.7639   0.3000  -3.5745  -2.1565   0.0000   0.0000     
+  1  1  3 12  -0.2500  69.1094   0.3000  -3.0983  -2.1565   0.0000   0.0000     
+  2  3 12  3  -0.4306   7.5000  -0.5000  -6.9948  -1.0000   0.0000   0.0000     
+  2  3 11  3   1.8627   9.7180  -1.0000  -7.2224  -1.0000   0.0000   0.0000     
+  1  3 11  3   2.5000  23.9443   1.0000  -3.2267  -1.0000   0.0000   0.0000     
+  1  1  3 11   0.9114  62.5039  -0.2389  -3.2976  -1.0000   0.0000   0.0000     
+  2  1  3 11   0.5000  35.0000   0.5000  -4.0000  -1.0000   0.0000   0.0000     
+  9    ! Nr of hydrogen bonds;at1;at2;at3;Rhb;Dehb;vhb1                         
+  3  2  3   2.1200  -3.5800   1.4500  19.5000                                   
+  3  2  4   2.0000  -6.0000   1.7976   3.0000                                   
+  4  2  3   1.2000  -2.0000   1.7976   3.0000                                   
+  4  2  4   1.2979  -6.0000   1.7976   3.0000                                   
+  3  2  5   1.5000  -2.0000   1.7976   3.0000                                   
+  4  2  5   1.5000  -2.0000   1.7976   3.0000                                   
+  5  2  3   1.5000  -2.0000   1.7976   3.0000                                   
+  5  2  4   1.5000  -2.0000   1.7976   3.0000                                   
+  5  2  5   1.5000  -2.0000   1.7976   3.0000                                   
diff --git a/data/benchmarks/water/ffield_acks2_300.water b/data/benchmarks/water/ffield_acks2_300.water
new file mode 100644
index 0000000000000000000000000000000000000000..d4d7b7b94777b6c456acb5bc1ad1e40dcfaa4da0
--- /dev/null
+++ b/data/benchmarks/water/ffield_acks2_300.water
@@ -0,0 +1,364 @@
+Reactive MD-force field: Water                                                  
+ 39       ! Number of general parameters                                        
+   50.0000 !Overcoordination parameter                                          
+    9.5469 !Overcoordination parameter                                          
+   26.5405 !Valency angle conjugation parameter                                 
+    1.7224 !Triple bond stabilisation parameter                                 
+    6.8702 !Triple bond stabilisation parameter                                 
+   60.4850 !C2-correction                                                       
+    1.0588 !Undercoordination parameter                                         
+    4.6000 !Triple bond stabilisation parameter                                 
+   12.1176 !Undercoordination parameter                                         
+   13.3056 !Undercoordination parameter                                         
+  -70.5044 !Triple bond stabilization energy                                    
+    0.0000 !Lower Taper-radius                                                  
+   10.0000 !Upper Taper-radius                                                  
+    2.8793 !Not used                                                            
+   33.8667 !Valency undercoordination                                           
+    6.0891 !Valency angle/lone pair parameter                                   
+    1.0563 !Valency angle                                                       
+    2.0384 !Valency angle parameter                                             
+    6.1431 !Not used                                                            
+    6.9290 !Double bond/angle parameter                                         
+    0.3989 !Double bond/angle parameter: overcoord                              
+    3.9954 !Double bond/angle parameter: overcoord                              
+   -2.4837 !Not used                                                            
+    5.7796 !Torsion/BO parameter                                                
+   10.0000 !Torsion overcoordination                                            
+    1.9487 !Torsion overcoordination                                            
+   -1.2327 !Conjugation 0 (not used)                                            
+    2.1645 !Conjugation                                                         
+    1.5591 !vdWaals shielding                                                   
+    0.1000 !Cutoff for bond order (*100)                                        
+    2.1365 !Valency angle conjugation parameter                                 
+    0.6991 !Overcoordination parameter                                          
+   50.0000 !Overcoordination parameter                                          
+    1.8512 !Valency/lone pair parameter                                         
+  548.6451 !Softness                                                            
+   20.0000 !Not used                                                            
+    5.0000 !Molecular energy (not used)                                         
+    0.0000 !Molecular energy (not used)                                         
+    2.6962 !Valency angle conjugation parameter                                 
+ 15    ! Nr of atoms; cov.r; valency;a.m;Rvdw;Evdw;gammaEEM;cov.r2;#            
+            alfa;gammavdW;valency;Eunder;Eover;chiEEM;etaEEM;n.u.               
+            cov r3;Elp;Heat inc.;n.u.;n.u.;n.u.;n.u.                            
+            ov/un;val1;n.u.;val3,vval4                                          
+ C    1.3817   4.0000  12.0000   1.8903   0.1838   0.9000   1.1341   4.0000     
+      9.7559   2.1346   4.0000  34.9350  79.5548   5.9666   7.0000   0.0000     
+      1.2114   0.0000 202.5551   8.9539  34.9289  13.5366   0.8563   0.0000     
+     -2.8983   2.5000   1.0564   4.0000   2.9663   0.0000   0.0000   0.0000     
+ H    0.8873   1.0000   1.0080   1.5420   0.0598   0.6883  -0.1000   1.0000     
+      8.1910  30.9706   1.0000   0.0000 121.1250   3.5768  10.5896   1.0000     
+     -0.1000   0.0000  61.6606   1.3986   2.1457   0.0003   3.4114   0.0000     
+    -15.7683   2.1488   1.0338   1.0000   2.8793   0.0000   0.0000   0.0000     
+ O    1.2450   2.0000  15.9990   2.3878   0.1023   1.0903   1.0548   6.0000     
+     10.5750  32.3923   4.0000  37.5000 116.0768   8.5000   7.5600   2.0000     
+      0.9049  -1.0100  59.0626   2.7162   3.2532   0.0021   0.9745   0.0000     
+     -3.6141   2.7025   1.0493   4.0000   2.9225   0.0000   5.4479   0.0000     
+ N    1.2333   3.0000  14.0000   1.9324   0.1376   0.8596   1.1748   5.0000     
+     10.0667   7.8431   4.0000  32.2482 100.0000   6.8418   6.3404   2.0000     
+      1.0433  13.7673 119.9837   2.1961   3.0696   2.7683   0.9745   0.0000     
+     -4.3875   2.6192   1.0183   4.0000   2.8793   0.0000   0.0000   0.0000     
+ S    1.9405   2.0000  32.0600   2.0677   0.2099   1.0336   1.5479   6.0000     
+      9.9575   4.9055   4.0000  52.9998 112.1416   6.5000   8.2545   2.0000     
+      1.4601   9.7177  71.1843   5.7487  23.2859  12.7147   0.9745   0.0000     
+    -11.0000   2.7466   1.0338   6.2998   2.8793   0.0000   0.0000   0.0000     
+ Si   2.0276   4.0000  28.0600   2.2042   0.1322   0.8218   1.5758   4.0000     
+     11.9413   2.0618   4.0000  11.8211 136.4845   1.8038   7.3852   0.0000     
+     -1.0000   0.0000 126.5182   6.4918   8.5961   0.2368   0.8563   0.0000     
+     -3.8112   3.1873   1.0338   6.2998   2.5791   0.0000   0.0000   0.0000     
+ Pt   1.9907   3.0000 195.0800   1.9980   0.2452   0.8218  -1.0000   3.0000     
+     12.8669   3.2118   3.0000   0.0000   0.0000   1.8038   7.3852   0.0000     
+     -1.0000   0.0000 142.6300   6.2293   5.2294   0.1542   0.8563   0.0000     
+     -6.7740   2.9867   1.0338   6.2998   2.5791   0.0000   0.0000   0.0000     
+ Zr   2.1000   4.0000  91.2240   2.1970   0.2542   0.8218  -1.0000   4.0000     
+     12.8545   3.5938   4.0000   0.0000   0.0000   1.8038   7.3852   0.0000     
+     -1.0000   0.0000 107.6300   6.2293   5.2294   0.1542   0.8563   0.0000     
+     -3.2224   2.9867   1.0338   6.2998   2.5791   0.0000   0.0000   0.0000     
+ Ni   1.8503   2.0000  58.6900   1.9219   0.1582   0.8218  -1.0000   2.0000     
+     12.1238   4.0351   2.0000   0.0000   0.0000   1.8038   7.3852   0.0000     
+     -1.0000   0.0000  95.6300   6.2293   5.2294   0.1542   0.8563   0.0000     
+     -3.2224   2.9867   1.0338   6.2998   2.5791   0.0000   0.0000   0.0000     
+ Au   1.8503   1.0000 196.9665   1.9219   0.1582   0.8218  -1.0000   1.0000     
+     12.1238   4.0351   1.0000   0.0000   0.0000   1.8038   7.3852   0.0000     
+     -1.0000   0.0000  72.6300   6.2293   5.2294   0.1542   0.8563   0.0000     
+     -3.2224   2.9867   1.0338   6.2998   2.5791   0.0000   0.0000   0.0000     
+ V    2.2657   3.0000  50.9415   1.7992   0.3005   0.6743   0.1000   5.0000     
+     12.3879   5.2243   3.0000   0.0000   0.0000  -0.3628   6.6023   0.0000     
+     -1.0000   0.0000 117.6300  23.1946   6.5795   0.0000   0.8563   0.0000     
+     -3.5389   1.5012   1.0338   3.0000   3.6411   0.0000   0.0000   0.0000     
+ Bi   2.1949   3.0000 208.9804   2.4429   0.1607   0.4960   0.0535   5.0000     
+     12.9571  35.5167   3.0000   0.0000   0.0000  -0.1926   6.4153   0.0000     
+     -1.0000   0.5785  52.6300   3.8978   0.9856   0.0314   0.8563   0.0000     
+     -2.5000   5.0597   1.0338   6.0000   2.5791   0.0000   0.0000   0.0000     
+ Ti   0.1000   4.0000  47.8800   2.0000   0.1659   0.6037   0.1000   4.0000     
+     13.2535   4.0063   4.0000  -5.0000   0.0000  -0.1864   5.9304   0.0000     
+     -1.0000   0.0000 129.6300  22.8461   1.8515   0.0064   0.8563   0.0000     
+     -3.4122   3.2711   1.0338   6.2998   2.2632   0.0000   0.0000   0.0000     
+ Mo   2.4710   5.6504  95.9400   1.8000   0.3285   1.0000   0.1000   6.0000     
+     13.0000  45.0000   4.0000   0.0000   0.0000   0.6062   6.1484   0.0000     
+      0.1000   0.0000 152.6300   3.7659   0.0689   2.9902   0.8563   0.0000     
+    -16.7660   3.1072   1.0338   8.0000   3.4590   0.0000   0.0000   0.0000     
+ X   -0.1000   2.0000   1.0080   2.0000   0.0000   1.0000  -0.1000   6.0000     
+     10.0000   2.5000   4.0000   0.0000   0.0000   8.5000   1.5000   0.0000     
+     -0.1000   0.0000  -2.3700   8.7410  13.3640   0.6690   0.9745   0.0000     
+    -11.0000   2.7466   1.0338   2.0000   2.8793   0.0000   0.0000   0.0000     
+ 40      ! Nr of bonds; Edis1;LPpen;n.u.;pbe1;pbo5;13corr;pbo6                  
+                         pbe2;pbo3;pbo4;n.u.;pbo1;pbo2;ovcorr                   
+  1  1 158.2004  99.1897  78.0000  -0.7738  -0.4550   1.0000  37.6117   0.4147  
+         0.4590  -0.1000   9.1628   1.0000  -0.0777   6.7268   1.0000   0.0000  
+  1  2 169.4760   0.0000   0.0000  -0.6083   0.0000   1.0000   6.0000   0.7652  
+         5.2290   1.0000   0.0000   1.0000  -0.0500   6.9136   0.0000   0.0000  
+  2  2 166.8880   0.0000   0.0000  -0.2191   0.0000   1.0000   6.0000   1.0000  
+         6.1152   1.0000   0.0000   1.0000  -0.0889   6.0000   0.0000   0.0000  
+  1  3 158.6946 107.4583  23.3136  -0.4240  -0.1743   1.0000  10.8209   1.0000  
+         0.5322  -0.3113   7.0000   1.0000  -0.1447   5.2450   0.0000   0.0000  
+  3  3 200.0000 230.7321  50.8293   0.2506  -0.1000   1.0000  29.7503   0.6051  
+         0.3451  -0.1055   9.0000   1.0000  -0.1225   5.5000   1.0000   0.0000  
+  1  4 134.1215 140.2179  79.9745   0.0163  -0.1428   1.0000  27.0617   0.2000  
+         0.1387  -0.3681   7.1611   1.0000  -0.1000   5.0825   1.0000   0.0000  
+  3  4 130.8596 169.4551  40.0000   0.3837  -0.1639   1.0000  35.0000   0.2000  
+         1.0000  -0.3579   7.0004   1.0000  -0.1193   6.8773   1.0000   0.0000  
+  4  4 157.9384  82.5526 152.5336   0.4010  -0.1034   1.0000  12.4261   0.5828  
+         0.1578  -0.1509  11.9186   1.0000  -0.0861   5.4271   1.0000   0.0000  
+  2  3 163.1043   0.0000   0.0000  -0.4155   0.0000   1.0000   6.0000   0.3607  
+         1.9380   1.0000   0.0000   0.0000  -0.0778   4.3082   0.0000   0.0000  
+  2  4 231.8173   0.0000   0.0000  -0.3364   0.0000   1.0000   6.0000   0.4402  
+         8.8910   1.0000   0.0000   1.0000  -0.0327   6.5754   0.0000   0.0000  
+  1  5 128.9942  74.5848  55.2528   0.1035  -0.5211   1.0000  18.9617   0.6000  
+         0.2949  -0.2398   8.1175   1.0000  -0.1029   5.6731   1.0000   0.0000  
+  2  5 151.5159   0.0000   0.0000  -0.4721   0.0000   1.0000   6.0000   0.6000  
+         9.4366   1.0000   0.0000   1.0000  -0.0290   7.0050   1.0000   0.0000  
+  3  5   0.0000   0.0000   0.0000   0.5563  -0.4038   1.0000  49.5611   0.6000  
+         0.4259  -0.4577  12.7569   1.0000  -0.1100   7.1145   1.0000   0.0000  
+  4  5   0.0000   0.0000   0.0000   0.4438  -0.2034   1.0000  40.3399   0.6000  
+         0.3296  -0.3153   9.1227   1.0000  -0.1805   5.6864   1.0000   0.0000  
+  5  5  96.1871  93.7006  68.6860   0.0955  -0.4781   1.0000  17.8574   0.6000  
+         0.2723  -0.2373   9.7875   1.0000  -0.0950   6.4757   1.0000   0.0000  
+  6  6 109.1904  70.8314  30.0000   0.2765  -0.3000   1.0000  16.0000   0.1583  
+         0.2804  -0.1994   8.1117   1.0000  -0.0675   8.2993   0.0000   0.0000  
+  2  6 137.1002   0.0000   0.0000  -0.1902   0.0000   1.0000   6.0000   0.4256  
+        17.7186   1.0000   0.0000   1.0000  -0.0377   6.4281   0.0000   0.0000  
+  3  6 191.1743  52.0733  43.3991  -0.2584  -0.3000   1.0000  36.0000   0.8764  
+         1.0248  -0.3658   4.2151   1.0000  -0.5004   4.2605   1.0000   0.0000  
+  4  6 185.4488  39.2832  43.3991  -0.1922  -0.3000   1.0000  36.0000   0.8217  
+         0.8538  -0.3887   4.4334   1.0000  -0.5241   4.4529   1.0000   0.0000  
+  7  7  90.1462   0.0000   0.0000   0.0004  -0.2000   0.0000  16.0000   0.3484  
+         1.0000  -0.2000  15.0000   1.0000  -0.1014   5.7631   0.0000   0.0000  
+  8  8  85.2900   0.0000   0.0000   0.0004  -0.2000   0.0000  16.0000   0.5438  
+         1.0000  -0.2000  15.0000   1.0000  -0.1001   5.5699   0.0000   0.0000  
+  9  9  73.6182   0.0000   0.0000   0.0004  -0.2000   0.0000  16.0000   0.3418  
+         1.0000  -0.2000  15.0000   1.0000  -0.1015   5.7850   0.0000   0.0000  
+ 10 10  73.6182   0.0000   0.0000   0.0004  -0.2000   0.0000  16.0000   0.3418  
+         1.0000  -0.2000  15.0000   1.0000  -0.1015   5.7850   0.0000   0.0000  
+ 11 11  36.2751   0.0000   0.0000   0.8059  -0.3000   0.0000  16.0000   0.1826  
+         0.3414  -0.3000  16.0000   1.0000  -0.0717   7.9108   0.0000   0.0000  
+  3 11 106.8008  67.5543   0.0000   0.0323  -0.3000   1.0000  36.0000   0.1000  
+         0.2670  -0.3402  16.0000   1.0000  -0.1761   4.6698   1.0000   0.0000  
+  2 11   0.0000   0.0000   0.0000  -0.2872  -0.3000   1.0000  36.0000   0.0082  
+         1.7973  -0.2500  20.0000   1.0000  -0.2578   6.5219   1.0000   0.0000  
+  1 11   0.0000   0.0000   0.0000  -0.2872  -0.3000   1.0000  36.0000   0.0082  
+         1.7973  -0.2500  20.0000   1.0000  -0.2578   6.5219   1.0000   0.0000  
+ 12 12  66.0677   0.0000   0.0000  -0.9557  -0.2000   0.0000  16.0000   0.2865  
+         0.5847  -0.2000  15.0000   1.0000  -0.0856   5.2857   0.0000   0.0000  
+  3 12 152.2407  57.6204   0.0000  -0.8033  -0.3000   1.0000  36.0000   0.0498  
+         1.8097  -0.3800  16.0000   1.0000  -0.2379   8.0000   1.0000   0.0000  
+  2 12  95.9209   0.0000   0.0000  -0.0153  -0.3000   1.0000  36.0000   0.0100  
+         1.0000  -0.2062   8.6647   1.0000  -0.1911   4.0000   1.0000   0.0000  
+  1 12  78.9091  40.6322   0.0000   0.0040  -0.3000   1.0000  36.0000   0.0384  
+         0.0904  -0.1209  12.3682   1.0000  -0.1613   4.3849   1.0000   0.0000  
+ 13 13  71.3016  10.0000   0.0000  -0.1571  -0.2000   0.0000  16.0000   0.3311  
+         0.1822  -0.2000  15.0000   1.0000  -0.1860   6.5172   0.0000   0.0000  
+  3 13 112.7130  29.8084   0.0000  -0.9010  -0.3000   1.0000  36.0000   0.5508  
+         0.1006  -0.2492  16.9476   1.0000  -0.1919   5.4797   1.0000   0.0000  
+  1 13   0.0000   0.0000   0.0000  -0.2872  -0.3000   1.0000  36.0000   0.0082  
+         1.7973  -0.2500  20.0000   1.0000  -0.2578   6.5219   1.0000   0.0000  
+  2 13   0.0000   0.0000   0.0000  -0.2872  -0.3000   1.0000  36.0000   0.0082  
+         1.7973  -0.2500  20.0000   1.0000  -0.2578   6.5219   1.0000   0.0000  
+  1 14   0.5356   0.9614   0.0000   0.3817  -0.3000   1.0000  36.0000   0.2142  
+         0.6116  -0.2579   6.1366   1.0000  -0.0913   6.6008   1.0000   0.0000  
+  2 14   0.0000   0.0000   0.0000  -0.2872  -0.3000   1.0000  36.0000   0.0082  
+         1.7973  -0.3027   4.6243   1.0000  -0.4578   3.5219   1.0000   0.0000  
+  3 14 112.7070  10.0000 135.5011   0.9277  -0.2354   1.0000  19.1731   1.2334  
+         0.9822  -0.1837   7.2216   1.0000  -0.1264   6.1257   1.0000   0.0000  
+ 14 14  44.6382   0.0000   0.0000   1.0000  -0.3000   0.0000  16.0000   0.2890  
+         0.3384  -0.3000  16.0000   1.0000  -0.1862   7.4588   0.0000   0.0000  
+ 12 14  50.0000   0.0000   0.0000   0.1000  -0.3000   0.0000  16.0000   0.3000  
+         1.0000  -0.3000  16.0000   1.0000  -0.2000   8.0000   0.0000   0.0000  
+ 20    ! Nr of off-diagonal terms; Ediss;Ro;gamma;rsigma;rpi;rpi2               
+  1  2   0.1239   1.4004   9.8467   1.1210  -1.0000  -1.0000                    
+  2  3   0.0299   1.3153  10.9102   0.9093  -1.0000  -1.0000                    
+  2  4   0.1059   1.8290   9.7818   0.9598  -1.0000  -1.0000                    
+  1  3   0.1156   1.8520   9.8317   1.2854   1.1352   1.0706                    
+  1  4   0.1447   1.8766   9.7990   1.3436   1.1885   1.1363                    
+  3  4   0.1048   2.0003  10.1220   1.3173   1.1096   1.0206                    
+  2  6   0.0470   1.6738  11.6877   1.1931  -1.0000  -1.0000                    
+  3  6   0.1263   1.8163  10.6833   1.6266   1.2052  -1.0000                    
+  1 11   0.1995   2.2133  13.0000   0.0102   1.4868  -1.0000                    
+  2 11   0.1319   1.5855  12.5457   0.0099   1.5065  -1.0000                    
+  3 11   0.0813   1.8649  10.8791   1.6498   1.6445  -1.0000                    
+  1 12   0.4235   1.7716  11.3664   1.8000   1.7212  -1.0000                    
+  2 12   0.0754   1.6033  12.4204   1.6896  -1.5000  -1.0000                    
+  3 12   0.1648   2.1260  11.2425   2.0692   1.6939  -1.0000                    
+  2 13   0.1340   1.8546  11.5784   1.0000  -1.0000  -1.0000                    
+  3 13   0.1280   1.8000  10.5743   1.7358   1.5296  -1.0000                    
+  1 13   0.1301   1.9382  11.1255   0.0100  -1.0000  -1.0000                    
+  1 14   0.1495   2.0794  12.2376   0.0100   1.4060  -1.0000                    
+  2 14   0.0795   1.6794  11.2376   0.0100   1.2060  -1.0000                    
+  3 14   0.2101   2.0342  10.4729   1.6019   1.4781   1.6548                    
+ 97    ! Nr of angles;at1;at2;at3;Thetao,o;ka;kb;pv1;pv2                        
+  1  1  1  59.0573  30.7029   0.7606   0.0000   0.7180   6.2933   1.1244        
+  1  1  2  65.7758  14.5234   6.2481   0.0000   0.5665   0.0000   1.6255        
+  2  1  2  70.2607  25.2202   3.7312   0.0000   0.0050   0.0000   2.7500        
+  1  2  2   0.0000   0.0000   6.0000   0.0000   0.0000   0.0000   1.0400        
+  1  2  1   0.0000   3.4110   7.7350   0.0000   0.0000   0.0000   1.0400        
+  2  2  2   0.0000  27.9213   5.8635   0.0000   0.0000   0.0000   1.0400        
+  1  1  3  49.6811   7.1713   4.3889   0.0000   0.7171  10.2661   1.0463        
+  3  1  3  77.7473  40.1718   2.9802 -25.3063   1.6170 -46.1315   2.2503        
+  1  1  4  66.1305  12.4661   7.0000   0.0000   3.0000  50.0000   1.1880        
+  3  1  4  73.9544  12.4661   7.0000   0.0000   3.0000   0.0000   1.1880        
+  4  1  4  64.1581  12.4661   7.0000   0.0000   3.0000   0.0000   1.1880        
+  2  1  3  65.0000  13.8815   5.0583   0.0000   0.4985   0.0000   1.4900        
+  2  1  4  74.2929  31.0883   2.6184   0.0000   0.0755   0.0000   1.0500        
+  1  2  4   0.0000   0.0019   6.3000   0.0000   0.0000   0.0000   1.0400        
+  1  3  1  73.5312  44.7275   0.7354   0.0000   3.0000   0.0000   1.0684        
+  1  3  3  79.4761  36.3701   1.8943   0.0000   0.7351  67.6777   3.0000        
+  1  3  4  82.4890  31.4554   0.9953   0.0000   1.6310   0.0000   1.0783        
+  3  3  3  80.7324  30.4554   0.9953   0.0000   1.6310  50.0000   1.0783        
+  3  3  4  84.3637  31.4554   0.9953   0.0000   1.6310   0.0000   1.0783        
+  4  3  4  89.7071  31.4554   0.9953   0.0000   1.6310   0.0000   1.1519        
+  1  3  2  70.1880  20.9562   0.3864   0.0000   0.0050   0.0000   1.6924        
+  2  3  3  75.6935  25.0000   2.0000   0.0000   1.0000   0.0000   1.1680        
+  2  3  4  75.6201  18.7919   0.9833   0.0000   0.1218   0.0000   1.0500        
+  2  3  2  77.3619   4.8342   7.1628   0.0000   2.9933   0.0000   1.5948        
+  1  4  1  66.0330  22.0295   1.4442   0.0000   1.6777   0.0000   1.0500        
+  1  4  3 103.3204  33.0381   0.5787   0.0000   1.6777   0.0000   1.0500        
+  1  4  4 104.1335   8.6043   1.6495   0.0000   1.6777   0.0000   1.0500        
+  3  4  3  74.1978  42.1786   1.7845 -18.0069   1.6777   0.0000   1.0500        
+  3  4  4  74.8600  43.7354   1.1572  -0.9193   1.6777   0.0000   1.0500        
+  4  4  4  75.0538  14.8267   5.2794   0.0000   1.6777   0.0000   1.0500        
+  1  4  2  69.1106  25.5067   1.1003   0.0000   0.0222   0.0000   1.0369        
+  2  4  3  81.3686  40.0712   2.2396   0.0000   0.0222   0.0000   1.0369        
+  2  4  4  83.0104  43.4766   1.5328   0.0000   0.0222   0.0000   1.0500        
+  2  4  2  70.8687  12.0168   5.0132   0.0000   0.0222   0.0000   1.1243        
+  1  2  3   0.0000  25.0000   3.0000   0.0000   1.0000   0.0000   1.0400        
+  1  2  4   0.0000   0.0019   6.0000   0.0000   0.0000   0.0000   1.0400        
+  1  2  5   0.0000   0.0019   6.0000   0.0000   0.0000   0.0000   1.0400        
+  3  2  3   0.0000   4.4124   2.5758   0.0000   0.0000   0.0000   2.9884        
+  3  2  4   0.0000   0.0019   6.0000   0.0000   0.0000   0.0000   1.0400        
+  4  2  4   0.0000   0.0019   6.0000   0.0000   0.0000   0.0000   1.0400        
+  2  2  3   0.0000  15.0000   2.2970   0.0000   0.0000   0.0000   1.3268        
+  2  2  4   0.0000   0.0019   6.0000   0.0000   0.0000   0.0000   1.0400        
+  1  1  5  74.9397  25.0560   1.8787   0.1463   0.0559   0.0000   1.0400        
+  1  5  1  86.9521  36.9951   2.0903   0.1463   0.0559   0.0000   1.0400        
+  2  1  5  74.9397  25.0560   1.8787   0.0000   0.0000   0.0000   1.0400        
+  1  5  2  86.1791  36.9951   2.0903   0.0000   0.0000   0.0000   1.0400        
+  1  5  5  85.3644  36.9951   2.0903   0.1463   0.0559   0.0000   1.0400        
+  2  5  2  93.1959  36.9951   2.0903   0.0000   0.0000   0.0000   1.0400        
+  2  5  5  84.3331  36.9951   2.0903   0.0000   0.0000   0.0000   1.0400        
+  6  6  6  69.3456  21.7361   1.4283   0.0000  -0.2101   0.0000   1.3241        
+  2  6  6  75.6168  21.5317   1.0435   0.0000   2.5179   0.0000   1.0400        
+  2  6  2  78.3939  20.9772   0.8630   0.0000   2.8421   0.0000   1.0400        
+  3  6  6  70.3016  15.4081   1.3267   0.0000   2.1459   0.0000   1.0400        
+  2  6  3  73.8232  16.6592   3.7425   0.0000   0.8613   0.0000   1.0400        
+  3  6  3  90.0344   7.7656   1.7264   0.0000   0.7689   0.0000   1.0400        
+  6  3  6  22.1715   3.6615   0.3160   0.0000   4.1125   0.0000   1.0400        
+  2  3  6  83.7634   5.6693   2.7780   0.0000   1.6982   0.0000   1.0400        
+  3  3  6  73.4663  25.0761   0.9143   0.0000   2.2466   0.0000   1.0400        
+  2  2  6   0.0000  47.1300   6.0000   0.0000   1.6371   0.0000   1.0400        
+  6  2  6   0.0000  31.5209   6.0000   0.0000   1.6371   0.0000   1.0400        
+  3  2  6   0.0000  31.0427   4.5625   0.0000   1.6371   0.0000   1.0400        
+  2  2  5   0.0000   0.0019   6.0000   0.0000   0.0000   0.0000   1.0400        
+  3 11  3  62.4906  31.5023   1.3328   0.0000   2.8731   0.0000   1.0794        
+ 11  3 11  31.0790  19.3435   0.4919   0.0000   2.9625   0.0000   3.0000        
+  3  3 11 100.0000  14.7642   7.0000   0.0000   1.0585   0.0000   1.1599        
+  1  3 11  60.7895  13.6681   0.7546   0.0000   2.1747   0.0000   2.9508        
+  2  3 11 100.0000   5.0000   1.4335   0.0000   1.2363   0.0000   5.0000        
+  3 12 12  23.8296   8.9089   7.0000   0.0000   1.0000   0.0000   2.8891        
+  3 12  3  87.0764  19.4489   2.5080   0.0000   2.6056   0.0000   3.0000        
+ 12  3 12  72.7369  13.7522   5.0243   0.0000   2.9700   0.0000   1.5506        
+  3  3 12  68.8771  10.5000   2.5500   0.0000   2.5729   0.0000   1.5892        
+  2  3 12  99.5836   5.4142   2.2105   0.0000   1.0513   0.0000   1.1000        
+  1  3 12  90.0000  12.1772   2.2055   0.0000   1.9064   0.0000   2.6056        
+  1  1 12  71.1708  32.6379   0.4516   0.0000   2.1609   0.0000   1.1000        
+  1 12  3  90.0000  45.0000   0.9335   0.0000   0.2140   0.0000   1.4846        
+  1 12  1  87.6204  45.0000   1.2740   0.0000   1.1519   0.0000   1.1000        
+  3  1 12  54.7020   3.2967   7.0000   0.0000   2.0408   0.0000   2.4032        
+  2 12  3  90.0000  28.2099   1.8036   0.0000   1.5461   0.0000   1.2304        
+  2 12  2  90.0000  36.3001   0.6409   0.0000   3.0000   0.0000   1.7755        
+  1 12  2  89.5835  45.0000   0.8465   0.0000   1.2118   0.0000   2.2282        
+  2  1 12  68.7714  22.9669   0.4631   0.0000   2.4269   0.0000   1.4680        
+  2  2 12   0.0000  30.2898   3.9181   0.0000   0.9914   0.0000   1.3121        
+  3  2 12   0.0000   1.0000   4.1706   0.0000   1.0100   0.0000   1.1000        
+  1  2 12   0.0000   1.0000   3.9722   0.0000   1.0075   0.0000   1.2984        
+  3 13  3  73.6321  10.6453   2.7693   0.0000   0.0500   0.0000   1.9906        
+ 13  3 13 100.0000   5.0270   5.0000   0.0000   1.2768   0.0000   2.0630        
+  3  3 13  52.3127  40.0000   1.1362   0.0000   1.5100   0.0000   1.1000        
+  3 13 13  66.6695   0.0036   3.2646   0.0000   0.0581   0.0000   1.3741        
+  2  3 13 100.0000   3.8927   8.0000   0.0000   2.0000   0.0000   1.1000        
+  1  3 13  96.6040   9.4537   8.0000   0.0000   0.3285   0.0000   4.0000        
+  3 14  3  79.6765  50.0000   1.0502  -0.0016   0.1000   0.0000   1.4583        
+ 14  3 14  20.2100  37.6165   0.6059   0.0000   0.1531   0.0000   2.0586        
+  3  3 14  38.5570  11.9307   0.9911   0.0000   0.8422   0.0000   1.0500        
+  3 14 14   5.8342   0.0724   0.1000   0.0000   0.5490   0.0000   1.7839        
+  2  3 14  81.8943   7.2820   2.1490   0.0000   0.6873   0.0000   3.2184        
+  1  3 14  75.5634   8.3289   1.0236   0.0000   2.0875   0.0000   1.0500        
+ 12  3 14  30.0000   5.0000   0.5000   0.0000   0.5000   0.0000   1.2500        
+ 47    ! Nr of torsions;at1;at2;at3;at4;;V1;V2;V3;V2(BO);vconj;n.u;n            
+  1  1  1  1  -0.2500  34.7453   0.0288  -6.3507  -1.6000   0.0000   0.0000     
+  1  1  1  2  -0.2500  29.2131   0.2945  -4.9581  -2.1802   0.0000   0.0000     
+  2  1  1  2  -0.2500  31.2081   0.4539  -4.8923  -2.2677   0.0000   0.0000     
+  1  1  1  3  -0.3495  22.2142  -0.2959  -2.5000  -1.9066   0.0000   0.0000     
+  2  1  1  3   0.0646  24.3195   0.6259  -3.9603  -1.0000   0.0000   0.0000     
+  3  1  1  3  -0.5456   5.5756   0.8433  -5.1924  -1.0180   0.0000   0.0000     
+  1  1  3  1   1.7555  27.9267   0.0072  -2.6533  -1.0000   0.0000   0.0000     
+  1  1  3  2  -1.4358  36.7830  -1.0000  -8.1821  -1.0000   0.0000   0.0000     
+  2  1  3  1  -1.3959  34.5053   0.7200  -2.5714  -2.1641   0.0000   0.0000     
+  2  1  3  2  -2.5000  70.0597   1.0000  -3.5539  -2.9929   0.0000   0.0000     
+  1  1  3  3   0.6852  11.2819  -0.4784  -2.5000  -2.1085   0.0000   0.0000     
+  2  1  3  3   0.1933  80.0000   1.0000  -4.0590  -3.0000   0.0000   0.0000     
+  3  1  3  1  -1.9889  76.4820  -0.1796  -3.8301  -3.0000   0.0000   0.0000     
+  3  1  3  2   0.2160  72.7707  -0.7087  -4.2100  -3.0000   0.0000   0.0000     
+  3  1  3  3  -2.5000  71.0772   0.2542  -3.1631  -3.0000   0.0000   0.0000     
+  1  3  3  1   2.5000  -0.6002   1.0000  -3.4297  -2.8858   0.0000   0.0000     
+  1  3  3  2  -2.5000  -3.3822   0.7004  -5.4467  -2.9586   0.0000   0.0000     
+  2  3  3  2   2.5000  -4.0000   0.9000  -2.5000  -1.0000   0.0000   0.0000     
+  1  3  3  3   1.2329  -4.0000   1.0000  -2.5000  -1.7479   0.0000   0.0000     
+  2  3  3  3   0.8302  -4.0000  -0.7763  -2.5000  -1.0000   0.0000   0.0000     
+  3  3  3  3  -2.5000  -4.0000   1.0000  -2.5000  -1.0000   0.0000   0.0000     
+  0  1  2  0   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000     
+  0  2  2  0   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000     
+  0  2  3  0   0.0000   0.1000   0.0200  -2.5415   0.0000   0.0000   0.0000     
+  0  1  1  0   0.0000  50.0000   0.3000  -4.0000  -2.0000   0.0000   0.0000     
+  0  3  3  0   0.5511  25.4150   1.1330  -5.1903  -1.0000   0.0000   0.0000     
+  0  1  4  0  -2.4242 128.1636   0.3739  -6.6098  -2.0000   0.0000   0.0000     
+  0  2  4  0   0.0000   0.1000   0.0200  -2.5415   0.0000   0.0000   0.0000     
+  0  3  4  0   1.4816  55.6641   0.0004  -7.0465  -2.7203   0.0000   0.0000     
+  0  4  4  0  -0.3244  27.7086   0.0039  -2.8272  -2.0000   0.0000   0.0000     
+  4  1  4  4  -5.5181   8.9706   0.0004  -6.1782  -2.0000   0.0000   0.0000     
+  0  1  5  0   3.3423  30.3435   0.0365  -2.7171   0.0000   0.0000   0.0000     
+  0  5  5  0  -0.0555 -42.7738   0.1515  -2.2056   0.0000   0.0000   0.0000     
+  0  2  5  0   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000     
+  0  6  6  0   0.0000   0.0000   0.1200  -2.4426   0.0000   0.0000   0.0000     
+  0  2  6  0   0.0000   0.0000   0.1200  -2.4847   0.0000   0.0000   0.0000     
+  0  3  6  0   0.0000   0.0000   0.1200  -2.4703   0.0000   0.0000   0.0000     
+  2  1  3 14   1.6297  56.8132   0.3398  -2.6912  -2.1000   0.0000   0.0000     
+  1  1  3 14  -0.0427  13.4096   0.9351  -6.5245  -2.1000   0.0000   0.0000     
+  2  3 14  3   2.5000  11.6208   1.0000  -9.0000  -1.0000   0.0000   0.0000     
+  2  1  3 12  -0.2500  45.7639   0.3000  -3.5745  -2.1565   0.0000   0.0000     
+  1  1  3 12  -0.2500  69.1094   0.3000  -3.0983  -2.1565   0.0000   0.0000     
+  2  3 12  3  -0.4306   7.5000  -0.5000  -6.9948  -1.0000   0.0000   0.0000     
+  2  3 11  3   1.8627   9.7180  -1.0000  -7.2224  -1.0000   0.0000   0.0000     
+  1  3 11  3   2.5000  23.9443   1.0000  -3.2267  -1.0000   0.0000   0.0000     
+  1  1  3 11   0.9114  62.5039  -0.2389  -3.2976  -1.0000   0.0000   0.0000     
+  2  1  3 11   0.5000  35.0000   0.5000  -4.0000  -1.0000   0.0000   0.0000     
+  9    ! Nr of hydrogen bonds;at1;at2;at3;Rhb;Dehb;vhb1                         
+  3  2  3   2.1845  -2.3549   3.0582  19.1627                                   
+  3  2  4   2.0000  -6.0000   1.7976   3.0000                                   
+  4  2  3   1.2000  -2.0000   1.7976   3.0000                                   
+  4  2  4   1.2979  -6.0000   1.7976   3.0000                                   
+  3  2  5   1.5000  -2.0000   1.7976   3.0000                                   
+  4  2  5   1.5000  -2.0000   1.7976   3.0000                                   
+  5  2  3   1.5000  -2.0000   1.7976   3.0000                                   
+  5  2  4   1.5000  -2.0000   1.7976   3.0000                                   
+  5  2  5   1.5000  -2.0000   1.7976   3.0000                                   
diff --git a/environ/param.gpu.water b/environ/param.gpu.water
index 3c38e57dd9a357e3d293dcedc2e0becf90f15e72..c7aaf52d7be169df93420f640d4abe7d6ce1aba0 100644
--- a/environ/param.gpu.water
+++ b/environ/param.gpu.water
@@ -59,3 +59,6 @@ freq_dipole_anal        1                       ! calculate electric dipole mome
 diffusion_coef          0                       ! 1: calculate diffusion coefficient of the system
 freq_diffusion_coef     1                       ! calculate diffusion coefficient at every 'this many' steps
 restrict_type           2                       ! -1: all types of atoms, 0 and up: only this type of atoms
+
+restart_format          1                       ! 0: restarts in ASCII  1: restarts in binary
+restart_freq            0                       ! 0: do not output any restart files. >0: output a restart file at every 'this many' steps
diff --git a/sPuReMD/src/charges.c b/sPuReMD/src/charges.c
index 54233396667390ad22f56163ceac981a645587c1..01279242f7a994bc9bd2f70b45a6742a45a2bde6 100644
--- a/sPuReMD/src/charges.c
+++ b/sPuReMD/src/charges.c
@@ -1420,7 +1420,7 @@ static void Extrapolate_Charges_QEq( const reax_system * const system,
 }
 
 
-static void Extrapolate_Charges_EEM( const reax_system * const system,
+static void Extrapolate_Charges_EE( const reax_system * const system,
         const control_params * const control,
         simulation_data * const data, static_storage * const workspace )
 {
@@ -1430,7 +1430,7 @@ static void Extrapolate_Charges_EEM( const reax_system * const system,
     /* extrapolation for s */
     //TODO: good candidate for vectorization, avoid moving data with head pointer and circular buffer
     #pragma omp parallel for schedule(static) \
-    default(none) private(i, s_tmp)
+        default(none) private(i, s_tmp)
     for ( i = 0; i < system->N_cm; ++i )
     {
         // no extrapolation
@@ -1466,6 +1466,7 @@ static void Compute_Preconditioner_QEq( const reax_system * const system,
         simulation_data * const data, static_storage * const workspace,
         const list * const far_nbrs )
 {
+    real time;
     sparse_matrix *Hptr;
 
     if ( control->cm_domain_sparsify_enabled == TRUE )
@@ -1477,6 +1478,22 @@ static void Compute_Preconditioner_QEq( const reax_system * const system,
         Hptr = workspace->H;
     }
 
+    time = Get_Time( );
+    if ( control->cm_solver_pre_app_type == TRI_SOLVE_GC_PA )
+    {
+        if ( control->cm_domain_sparsify_enabled == TRUE )
+        {
+            Hptr = setup_graph_coloring( workspace->H_sp );
+        }
+        else
+        {
+            Hptr = setup_graph_coloring( workspace->H );
+        }
+
+        Sort_Matrix_Rows( Hptr );
+    }
+    data->timing.cm_sort_mat_rows += Get_Timing_Info( time );
+
 #if defined(TEST_MAT)
     Hptr = create_test_mat( );
 #endif
@@ -1527,53 +1544,300 @@ static void Compute_Preconditioner_QEq( const reax_system * const system,
 
 #if defined(DEBUG_FOCUS)
         sprintf( fname, "%s.L%d.out", control->sim_name, data->step );
-        Print_Sparse_Matrix2( workspace->L, fname );
+        Print_Sparse_Matrix2( workspace->L, fname, NULL );
         sprintf( fname, "%s.U%d.out", control->sim_name, data->step );
-        Print_Sparse_Matrix2( workspace->U, fname );
-
-        fprintf( stderr, "icholt-" );
-        //sprintf( fname, "%s.L%d.out", control->sim_name, data->step );
-        //Print_Sparse_Matrix2( workspace->L, fname );
-        //Print_Sparse_Matrix( U );
+        Print_Sparse_Matrix2( workspace->U, fname, NULL );
 #endif
     }
 #endif
 }
 
 
-/* Compute preconditioner for EEM
+/* Compute preconditioner for EE
  */
-static void Compute_Preconditioner_EEM( const reax_system * const system,
+//static void Compute_Preconditioner_EE( const reax_system * const system,
+//        const control_params * const control,
+//        simulation_data * const data, static_storage * const workspace,
+//        const list * const far_nbrs )
+//{
+//    int i, top;
+//    static real * ones = NULL, * x = NULL, * y = NULL;
+//    sparse_matrix *Hptr;
+//
+//    Hptr = workspace->H_EE;
+//
+//#if defined(TEST_MAT)
+//    Hptr = create_test_mat( );
+//#endif
+//
+//    if ( ones == NULL )
+//    {
+//        if ( ( ones = (real*) malloc( system->N * sizeof(real)) ) == NULL ||
+//            ( x = (real*) malloc( system->N * sizeof(real)) ) == NULL ||
+//            ( y = (real*) malloc( system->N * sizeof(real)) ) == NULL )
+//        {
+//            fprintf( stderr, "Not enough space for preconditioner computation. Terminating...\n" );
+//            exit( INSUFFICIENT_MEMORY );
+//        }
+//
+//        for ( i = 0; i < system->N; ++i )
+//        {
+//            ones[i] = 1.0;
+//        }
+//    }
+//
+//    switch ( control->cm_solver_pre_comp_type )
+//    {
+//    case DIAG_PC:
+//        data->timing.cm_solver_pre_comp +=
+//            diag_pre_comp( Hptr, workspace->Hdia_inv );
+//        break;
+//
+//    case ICHOLT_PC:
+//        data->timing.cm_solver_pre_comp +=
+//            ICHOLT( Hptr, workspace->droptol, workspace->L_EE, workspace->U_EE );
+//        break;
+//
+//    case ILU_PAR_PC:
+//        data->timing.cm_solver_pre_comp +=
+//            ILU_PAR( Hptr, control->cm_solver_pre_comp_sweeps, workspace->L_EE, workspace->U_EE );
+//        break;
+//
+//    case ILUT_PAR_PC:
+//        data->timing.cm_solver_pre_comp +=
+//            ILUT_PAR( Hptr, workspace->droptol, control->cm_solver_pre_comp_sweeps,
+//                    workspace->L_EE, workspace->U_EE );
+//        break;
+//
+//    case ILU_SUPERLU_MT_PC:
+//#if defined(HAVE_SUPERLU_MT)
+//        data->timing.cm_solver_pre_comp +=
+//            SuperLU_Factorize( Hptr, workspace->L_EE, workspace->U_EE );
+//#else
+//        fprintf( stderr, "SuperLU MT support disabled. Re-compile before enabling. Terminating...\n" );
+//        exit( INVALID_INPUT );
+//#endif
+//        break;
+//
+//    default:
+//        fprintf( stderr, "Unrecognized preconditioner computation method. Terminating...\n" );
+//        exit( INVALID_INPUT );
+//        break;
+//    }
+//
+//    if ( control->cm_solver_pre_comp_type != DIAG_PC )
+//    {
+//        switch ( control->cm_solver_pre_app_type )
+//        {
+//            case TRI_SOLVE_PA:
+//                tri_solve( workspace->L_EE, ones, x, workspace->L_EE->n, LOWER );
+//                Transpose_I( workspace->U_EE );
+//                tri_solve( workspace->U_EE, ones, y, workspace->U_EE->n, LOWER );
+//                Transpose_I( workspace->U_EE );
+//
+//                memcpy( workspace->L->start, workspace->L_EE->start, sizeof(unsigned int) * (system->N + 1) );
+//                memcpy( workspace->L->j, workspace->L_EE->j, sizeof(unsigned int) * workspace->L_EE->start[workspace->L_EE->n] );
+//                memcpy( workspace->L->val, workspace->L_EE->val, sizeof(real) * workspace->L_EE->start[workspace->L_EE->n] );
+//
+//                top = workspace->L->start[system->N];
+//                for ( i = 0; i < system->N; ++i )
+//                {
+//                    workspace->L->j[top] = i;
+//                    workspace->L->val[top] = x[i];
+//                    ++top;
+//                }
+//
+//                workspace->L->j[top] = system->N_cm - 1;
+//                workspace->L->val[top] = 1.0;
+//                ++top;
+//
+//                workspace->L->start[system->N_cm] = top;
+//
+//                top = 0;
+//                for ( i = 0; i < system->N; ++i )
+//                {
+//                    workspace->U->start[i] = top;
+//                    memcpy( workspace->U->j + top, workspace->U_EE->j + workspace->U_EE->start[i],
+//                            sizeof(unsigned int) * (workspace->U_EE->start[i + 1] - workspace->U_EE->start[i]) );
+//                    memcpy( workspace->U->val + top, workspace->U_EE->val + workspace->U_EE->start[i],
+//                            sizeof(real) * (workspace->U_EE->start[i + 1] - workspace->U_EE->start[i]) );
+//                    top += (workspace->U_EE->start[i + 1] - workspace->U_EE->start[i]);
+//
+//                    workspace->U->j[top] = system->N_cm - 1;
+//                    workspace->U->val[top] = y[i];
+//                    ++top;
+//                }
+//
+//                workspace->U->start[system->N_cm - 1] = top;
+//
+//                workspace->U->j[top] = system->N_cm - 1;
+//                workspace->U->val[top] = -Dot( x, y, system->N );
+//                ++top;
+//
+//                workspace->U->start[system->N_cm] = top;
+//                break;
+//
+//            case TRI_SOLVE_LEVEL_SCHED_PA:
+//                tri_solve_level_sched( workspace->L_EE, ones, x, workspace->L_EE->n, LOWER, TRUE );
+//                Transpose_I( workspace->U_EE );
+//                tri_solve_level_sched( workspace->U_EE, ones, y, workspace->U_EE->n, LOWER, TRUE );
+//                Transpose_I( workspace->U_EE );
+//
+//                memcpy( workspace->L->start, workspace->L_EE->start, sizeof(unsigned int) * (system->N + 1) );
+//                memcpy( workspace->L->j, workspace->L_EE->j, sizeof(unsigned int) * workspace->L_EE->start[workspace->L_EE->n] );
+//                memcpy( workspace->L->val, workspace->L_EE->val, sizeof(real) * workspace->L_EE->start[workspace->L_EE->n] );
+//
+//                top = workspace->L->start[system->N];
+//                for ( i = 0; i < system->N; ++i )
+//                {
+//                    workspace->L->j[top] = i;
+//                    workspace->L->val[top] = x[i];
+//                    ++top;
+//                }
+//
+//                workspace->L->j[top] = system->N_cm - 1;
+//                workspace->L->val[top] = 1.0;
+//                ++top;
+//
+//                workspace->L->start[system->N_cm] = top;
+//
+//                top = 0;
+//                for ( i = 0; i < system->N; ++i )
+//                {
+//                    workspace->U->start[i] = top;
+//                    memcpy( workspace->U->j + top, workspace->U_EE->j + workspace->U_EE->start[i],
+//                            sizeof(unsigned int) * (workspace->U_EE->start[i + 1] - workspace->U_EE->start[i]) );
+//                    memcpy( workspace->U->val + top, workspace->U_EE->val + workspace->U_EE->start[i],
+//                            sizeof(real) * (workspace->U_EE->start[i + 1] - workspace->U_EE->start[i]) );
+//                    top += (workspace->U_EE->start[i + 1] - workspace->U_EE->start[i]);
+//
+//                    workspace->U->j[top] = system->N_cm - 1;
+//                    workspace->U->val[top] = y[i];
+//                    ++top;
+//                }
+//
+//                workspace->U->start[system->N_cm - 1] = top;
+//
+//                workspace->U->j[top] = system->N_cm - 1;
+//                workspace->U->val[top] = -Dot( x, y, system->N );
+//                ++top;
+//
+//                workspace->U->start[system->N_cm] = top;
+//                break;
+//
+//            //TODO: add Jacobi iter, etc.?
+//            default:
+//                tri_solve( workspace->L_EE, ones, x, workspace->L_EE->n, LOWER );
+//                Transpose_I( workspace->U_EE );
+//                tri_solve( workspace->U_EE, ones, y, workspace->U_EE->n, LOWER );
+//                Transpose_I( workspace->U_EE );
+//
+//                memcpy( workspace->L->start, workspace->L_EE->start, sizeof(unsigned int) * (system->N + 1) );
+//                memcpy( workspace->L->j, workspace->L_EE->j, sizeof(unsigned int) * workspace->L_EE->start[workspace->L_EE->n] );
+//                memcpy( workspace->L->val, workspace->L_EE->val, sizeof(real) * workspace->L_EE->start[workspace->L_EE->n] );
+//
+//                top = workspace->L->start[system->N];
+//                for ( i = 0; i < system->N; ++i )
+//                {
+//                    workspace->L->j[top] = i;
+//                    workspace->L->val[top] = x[i];
+//                    ++top;
+//                }
+//
+//                workspace->L->j[top] = system->N_cm - 1;
+//                workspace->L->val[top] = 1.0;
+//                ++top;
+//
+//                workspace->L->start[system->N_cm] = top;
+//
+//                top = 0;
+//                for ( i = 0; i < system->N; ++i )
+//                {
+//                    workspace->U->start[i] = top;
+//                    memcpy( workspace->U->j + top, workspace->U_EE->j + workspace->U_EE->start[i],
+//                            sizeof(unsigned int) * (workspace->U_EE->start[i + 1] - workspace->U_EE->start[i]) );
+//                    memcpy( workspace->U->val + top, workspace->U_EE->val + workspace->U_EE->start[i],
+//                            sizeof(real) * (workspace->U_EE->start[i + 1] - workspace->U_EE->start[i]) );
+//                    top += (workspace->U_EE->start[i + 1] - workspace->U_EE->start[i]);
+//
+//                    workspace->U->j[top] = system->N_cm - 1;
+//                    workspace->U->val[top] = y[i];
+//                    ++top;
+//                }
+//
+//                workspace->U->start[system->N_cm - 1] = top;
+//
+//                workspace->U->j[top] = system->N_cm - 1;
+//                workspace->U->val[top] = -Dot( x, y, system->N );
+//                ++top;
+//
+//                workspace->U->start[system->N_cm] = top;
+//                break;
+//        }
+//    }
+//
+//#if defined(DEBUG)
+//    if ( control->cm_solver_pre_comp_type != DIAG_PC )
+//    {
+//        fprintf( stderr, "condest = %f\n", condest(workspace->L) );
+//
+//#if defined(DEBUG_FOCUS)
+//        sprintf( fname, "%s.L%d.out", control->sim_name, data->step );
+//        Print_Sparse_Matrix2( workspace->L, fname, NULL );
+//        sprintf( fname, "%s.U%d.out", control->sim_name, data->step );
+//        Print_Sparse_Matrix2( workspace->U, fname, NULL );
+//
+//        fprintf( stderr, "icholt-" );
+//        sprintf( fname, "%s.L%d.out", control->sim_name, data->step );
+//        Print_Sparse_Matrix2( workspace->L, fname, NULL );
+//        Print_Sparse_Matrix( U );
+//#endif
+//    }
+//#endif
+//}
+
+
+/* Compute preconditioner for EE
+ */
+static void Compute_Preconditioner_EE( const reax_system * const system,
         const control_params * const control,
         simulation_data * const data, static_storage * const workspace,
         const list * const far_nbrs )
 {
-    int i, top;
-    static real * ones = NULL, * x = NULL, * y = NULL;
+    real time;
     sparse_matrix *Hptr;
 
-    Hptr = workspace->H_EEM;
-
-#if defined(TEST_MAT)
-    Hptr = create_test_mat( );
-#endif
+    if ( control->cm_domain_sparsify_enabled == TRUE )
+    {
+        Hptr = workspace->H_sp;
+    }
+    else
+    {
+        Hptr = workspace->H;
+    }
 
-    if ( ones == NULL )
+    time = Get_Time( );
+    if ( control->cm_solver_pre_app_type == TRI_SOLVE_GC_PA )
     {
-        if ( ( ones = (real*) malloc( system->N * sizeof(real)) ) == NULL ||
-            ( x = (real*) malloc( system->N * sizeof(real)) ) == NULL ||
-            ( y = (real*) malloc( system->N * sizeof(real)) ) == NULL )
+        if ( control->cm_domain_sparsify_enabled == TRUE )
         {
-            fprintf( stderr, "Not enough space for preconditioner computation. Terminating...\n" );
-            exit( INSUFFICIENT_MEMORY );
+            Hptr = setup_graph_coloring( workspace->H_sp );
         }
-
-        for ( i = 0; i < system->N; ++i )
+        else
         {
-            ones[i] = 1.0;
+            Hptr = setup_graph_coloring( workspace->H );
         }
+
+        Sort_Matrix_Rows( Hptr );
     }
+    data->timing.cm_sort_mat_rows += Get_Timing_Info( time );
 
+#if defined(TEST_MAT)
+    Hptr = create_test_mat( );
+#endif
+
+    Hptr->val[Hptr->start[system->N + 1] - 1] = 1.0;
+    
     switch ( control->cm_solver_pre_comp_type )
     {
     case DIAG_PC:
@@ -1583,24 +1847,24 @@ static void Compute_Preconditioner_EEM( const reax_system * const system,
 
     case ICHOLT_PC:
         data->timing.cm_solver_pre_comp +=
-            ICHOLT( Hptr, workspace->droptol, workspace->L_EEM, workspace->U_EEM );
+            ICHOLT( Hptr, workspace->droptol, workspace->L, workspace->U );
         break;
 
     case ILU_PAR_PC:
         data->timing.cm_solver_pre_comp +=
-            ILU_PAR( Hptr, control->cm_solver_pre_comp_sweeps, workspace->L_EEM, workspace->U_EEM );
+            ILU_PAR( Hptr, control->cm_solver_pre_comp_sweeps, workspace->L, workspace->U );
         break;
 
     case ILUT_PAR_PC:
         data->timing.cm_solver_pre_comp +=
             ILUT_PAR( Hptr, workspace->droptol, control->cm_solver_pre_comp_sweeps,
-                    workspace->L_EEM, workspace->U_EEM );
+                    workspace->L, workspace->U );
         break;
 
     case ILU_SUPERLU_MT_PC:
 #if defined(HAVE_SUPERLU_MT)
         data->timing.cm_solver_pre_comp +=
-            SuperLU_Factorize( Hptr, workspace->L_EEM, workspace->U_EEM );
+            SuperLU_Factorize( Hptr, workspace->L, workspace->U );
 #else
         fprintf( stderr, "SuperLU MT support disabled. Re-compile before enabling. Terminating...\n" );
         exit( INVALID_INPUT );
@@ -1613,172 +1877,18 @@ static void Compute_Preconditioner_EEM( const reax_system * const system,
         break;
     }
 
-    if ( control->cm_solver_pre_comp_type != DIAG_PC )
-    {
-        switch ( control->cm_solver_pre_app_type )
-        {
-            case TRI_SOLVE_PA:
-                tri_solve( workspace->L_EEM, ones, x, workspace->L_EEM->n, LOWER );
-                Transpose_I( workspace->U_EEM );
-                tri_solve( workspace->U_EEM, ones, y, workspace->U_EEM->n, LOWER );
-                Transpose_I( workspace->U_EEM );
-
-                memcpy( workspace->L->start, workspace->L_EEM->start, sizeof(unsigned int) * (system->N + 1) );
-                memcpy( workspace->L->j, workspace->L_EEM->j, sizeof(unsigned int) * workspace->L_EEM->start[workspace->L_EEM->n] );
-                memcpy( workspace->L->val, workspace->L_EEM->val, sizeof(real) * workspace->L_EEM->start[workspace->L_EEM->n] );
-
-                top = workspace->L->start[system->N];
-                for ( i = 0; i < system->N; ++i )
-                {
-                    workspace->L->j[top] = i;
-                    workspace->L->val[top] = x[i];
-                    ++top;
-                }
-
-                workspace->L->j[top] = system->N_cm - 1;
-                workspace->L->val[top] = 1.0;
-                ++top;
-
-                workspace->L->start[system->N_cm] = top;
-
-                top = 0;
-                for ( i = 0; i < system->N; ++i )
-                {
-                    workspace->U->start[i] = top;
-                    memcpy( workspace->U->j + top, workspace->U_EEM->j + workspace->U_EEM->start[i],
-                            sizeof(unsigned int) * (workspace->U_EEM->start[i + 1] - workspace->U_EEM->start[i]) );
-                    memcpy( workspace->U->val + top, workspace->U_EEM->val + workspace->U_EEM->start[i],
-                            sizeof(real) * (workspace->U_EEM->start[i + 1] - workspace->U_EEM->start[i]) );
-                    top += (workspace->U_EEM->start[i + 1] - workspace->U_EEM->start[i]);
-
-                    workspace->U->j[top] = system->N_cm - 1;
-                    workspace->U->val[top] = y[i];
-                    ++top;
-                }
-
-                workspace->U->start[system->N_cm - 1] = top;
-
-                workspace->U->j[top] = system->N_cm - 1;
-                workspace->U->val[top] = -Dot( x, y, system->N );
-                ++top;
-
-                workspace->U->start[system->N_cm] = top;
-                break;
-
-            case TRI_SOLVE_LEVEL_SCHED_PA:
-                tri_solve_level_sched( workspace->L_EEM, ones, x, workspace->L_EEM->n, LOWER, TRUE );
-                Transpose_I( workspace->U_EEM );
-                tri_solve_level_sched( workspace->U_EEM, ones, y, workspace->U_EEM->n, LOWER, TRUE );
-                Transpose_I( workspace->U_EEM );
-
-                memcpy( workspace->L->start, workspace->L_EEM->start, sizeof(unsigned int) * (system->N + 1) );
-                memcpy( workspace->L->j, workspace->L_EEM->j, sizeof(unsigned int) * workspace->L_EEM->start[workspace->L_EEM->n] );
-                memcpy( workspace->L->val, workspace->L_EEM->val, sizeof(real) * workspace->L_EEM->start[workspace->L_EEM->n] );
-
-                top = workspace->L->start[system->N];
-                for ( i = 0; i < system->N; ++i )
-                {
-                    workspace->L->j[top] = i;
-                    workspace->L->val[top] = x[i];
-                    ++top;
-                }
-
-                workspace->L->j[top] = system->N_cm - 1;
-                workspace->L->val[top] = 1.0;
-                ++top;
-
-                workspace->L->start[system->N_cm] = top;
-
-                top = 0;
-                for ( i = 0; i < system->N; ++i )
-                {
-                    workspace->U->start[i] = top;
-                    memcpy( workspace->U->j + top, workspace->U_EEM->j + workspace->U_EEM->start[i],
-                            sizeof(unsigned int) * (workspace->U_EEM->start[i + 1] - workspace->U_EEM->start[i]) );
-                    memcpy( workspace->U->val + top, workspace->U_EEM->val + workspace->U_EEM->start[i],
-                            sizeof(real) * (workspace->U_EEM->start[i + 1] - workspace->U_EEM->start[i]) );
-                    top += (workspace->U_EEM->start[i + 1] - workspace->U_EEM->start[i]);
-
-                    workspace->U->j[top] = system->N_cm - 1;
-                    workspace->U->val[top] = y[i];
-                    ++top;
-                }
-
-                workspace->U->start[system->N_cm - 1] = top;
-
-                workspace->U->j[top] = system->N_cm - 1;
-                workspace->U->val[top] = -Dot( x, y, system->N );
-                ++top;
-
-                workspace->U->start[system->N_cm] = top;
-                break;
-
-            //TODO: add Jacobi iter, etc.?
-            default:
-                tri_solve( workspace->L_EEM, ones, x, workspace->L_EEM->n, LOWER );
-                Transpose_I( workspace->U_EEM );
-                tri_solve( workspace->U_EEM, ones, y, workspace->U_EEM->n, LOWER );
-                Transpose_I( workspace->U_EEM );
-
-                memcpy( workspace->L->start, workspace->L_EEM->start, sizeof(unsigned int) * (system->N + 1) );
-                memcpy( workspace->L->j, workspace->L_EEM->j, sizeof(unsigned int) * workspace->L_EEM->start[workspace->L_EEM->n] );
-                memcpy( workspace->L->val, workspace->L_EEM->val, sizeof(real) * workspace->L_EEM->start[workspace->L_EEM->n] );
-
-                top = workspace->L->start[system->N];
-                for ( i = 0; i < system->N; ++i )
-                {
-                    workspace->L->j[top] = i;
-                    workspace->L->val[top] = x[i];
-                    ++top;
-                }
-
-                workspace->L->j[top] = system->N_cm - 1;
-                workspace->L->val[top] = 1.0;
-                ++top;
-
-                workspace->L->start[system->N_cm] = top;
-
-                top = 0;
-                for ( i = 0; i < system->N; ++i )
-                {
-                    workspace->U->start[i] = top;
-                    memcpy( workspace->U->j + top, workspace->U_EEM->j + workspace->U_EEM->start[i],
-                            sizeof(unsigned int) * (workspace->U_EEM->start[i + 1] - workspace->U_EEM->start[i]) );
-                    memcpy( workspace->U->val + top, workspace->U_EEM->val + workspace->U_EEM->start[i],
-                            sizeof(real) * (workspace->U_EEM->start[i + 1] - workspace->U_EEM->start[i]) );
-                    top += (workspace->U_EEM->start[i + 1] - workspace->U_EEM->start[i]);
-
-                    workspace->U->j[top] = system->N_cm - 1;
-                    workspace->U->val[top] = y[i];
-                    ++top;
-                }
-
-                workspace->U->start[system->N_cm - 1] = top;
-
-                workspace->U->j[top] = system->N_cm - 1;
-                workspace->U->val[top] = -Dot( x, y, system->N );
-                ++top;
-
-                workspace->U->start[system->N_cm] = top;
-                break;
-        }
-    }
+    Hptr->val[Hptr->start[system->N + 1] - 1] = 0.0;
 
 #if defined(DEBUG)
     if ( control->cm_solver_pre_comp_type != DIAG_PC )
     {
-        fprintf( stderr, "condest = %f\n", condest(workspace->L) );
+        fprintf( stderr, "condest = %f\n", condest(workspace->L, workspace->U) );
 
 #if defined(DEBUG_FOCUS)
         sprintf( fname, "%s.L%d.out", control->sim_name, data->step );
-        Print_Sparse_Matrix2( workspace->L, fname );
+        Print_Sparse_Matrix2( workspace->L, fname, NULL );
         sprintf( fname, "%s.U%d.out", control->sim_name, data->step );
-        Print_Sparse_Matrix2( workspace->U, fname );
-
-        fprintf( stderr, "icholt-" );
-        //sprintf( fname, "%s.L%d.out", control->sim_name, data->step );
-        //Print_Sparse_Matrix2( workspace->L, fname );
-        //Print_Sparse_Matrix( U );
+        Print_Sparse_Matrix2( workspace->U, fname, NULL );
 #endif
     }
 #endif
@@ -1792,6 +1902,7 @@ static void Compute_Preconditioner_ACKS2( const reax_system * const system,
         simulation_data * const data, static_storage * const workspace,
         const list * const far_nbrs )
 {
+    real time;
     sparse_matrix *Hptr;
 
     if ( control->cm_domain_sparsify_enabled == TRUE )
@@ -1803,17 +1914,28 @@ static void Compute_Preconditioner_ACKS2( const reax_system * const system,
         Hptr = workspace->H;
     }
 
+    time = Get_Time( );
+    if ( control->cm_solver_pre_app_type == TRI_SOLVE_GC_PA )
+    {
+        if ( control->cm_domain_sparsify_enabled == TRUE )
+        {
+            Hptr = setup_graph_coloring( workspace->H_sp );
+        }
+        else
+        {
+            Hptr = setup_graph_coloring( workspace->H );
+        }
+
+        Sort_Matrix_Rows( Hptr );
+    }
+    data->timing.cm_sort_mat_rows += Get_Timing_Info( time );
+
 #if defined(TEST_MAT)
     Hptr = create_test_mat( );
 #endif
 
-//    fprintf( stderr, "(%4d, %4d): %f\n", system->N, Hptr->j[Hptr->start[system->N + 1] - 1], Hptr->val[Hptr->start[system->N + 1] - 1] );
-//    fprintf( stderr, "(%4d, %4d): %f\n", system->N_cm - 1, Hptr->j[Hptr->start[system->N_cm] - 1], Hptr->val[Hptr->start[system->N_cm] - 1] );
     Hptr->val[Hptr->start[system->N + 1] - 1] = 1.0;
     Hptr->val[Hptr->start[system->N_cm] - 1] = 1.0;
-//    fprintf( stderr, "(%4d, %4d): %f\n", system->N, Hptr->j[Hptr->start[system->N + 1] - 1], Hptr->val[Hptr->start[system->N + 1] - 1] );
-//    fprintf( stderr, "(%4d, %4d): %f\n", system->N_cm - 1, Hptr->j[Hptr->start[system->N_cm] - 1], Hptr->val[Hptr->start[system->N_cm] - 1] );
-//    fflush( stderr );
     
     switch ( control->cm_solver_pre_comp_type )
     {
@@ -1864,14 +1986,9 @@ static void Compute_Preconditioner_ACKS2( const reax_system * const system,
 
 #if defined(DEBUG_FOCUS)
         sprintf( fname, "%s.L%d.out", control->sim_name, data->step );
-        Print_Sparse_Matrix2( workspace->L, fname );
+        Print_Sparse_Matrix2( workspace->L, fname, NULL );
         sprintf( fname, "%s.U%d.out", control->sim_name, data->step );
-        Print_Sparse_Matrix2( workspace->U, fname );
-
-        fprintf( stderr, "icholt-" );
-        //sprintf( fname, "%s.L%d.out", control->sim_name, data->step );
-        //Print_Sparse_Matrix2( workspace->L, fname );
-        //Print_Sparse_Matrix( U );
+        Print_Sparse_Matrix2( workspace->U, fname, NULL );
 #endif
     }
 #endif
@@ -1905,20 +2022,6 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
     {
         Sort_Matrix_Rows( workspace->H_sp );
     }
-
-    if ( control->cm_solver_pre_app_type == TRI_SOLVE_GC_PA )
-    {
-        if ( control->cm_domain_sparsify_enabled == TRUE )
-        {
-            Hptr = setup_graph_coloring( workspace->H_sp );
-        }
-        else
-        {
-            Hptr = setup_graph_coloring( workspace->H );
-        }
-
-        Sort_Matrix_Rows( Hptr );
-    }
     data->timing.cm_sort_mat_rows += Get_Timing_Info( time );
 
 #if defined(DEBUG)
@@ -2034,9 +2137,9 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
 }
 
 
-/* Setup routines before computing the preconditioner for EEM
+/* Setup routines before computing the preconditioner for EE
  */
-static void Setup_Preconditioner_EEM( const reax_system * const system,
+static void Setup_Preconditioner_EE( const reax_system * const system,
         const control_params * const control,
         simulation_data * const data, static_storage * const workspace,
         const list * const far_nbrs )
@@ -2061,45 +2164,14 @@ static void Setup_Preconditioner_EEM( const reax_system * const system,
     {
         Sort_Matrix_Rows( workspace->H_sp );
     }
-
-    if ( control->cm_solver_pre_app_type == TRI_SOLVE_GC_PA )
-    {
-        if ( control->cm_domain_sparsify_enabled == TRUE )
-        {
-            Hptr = setup_graph_coloring( workspace->H_sp );
-        }
-        else
-        {
-            Hptr = setup_graph_coloring( workspace->H );
-        }
-
-        Sort_Matrix_Rows( Hptr );
-    }
     data->timing.cm_sort_mat_rows += Get_Timing_Info( time );
 
+    Hptr->val[Hptr->start[system->N + 1] - 1] = 1.0;
+
 #if defined(DEBUG)
     fprintf( stderr, "H matrix sorted\n" );
 #endif
 
-    if ( workspace->H_EEM == NULL )
-    {
-        if ( Allocate_Matrix( &(workspace->H_EEM), system->N, workspace->H->m ) == FAILURE )
-        {
-            fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
-            exit( INSUFFICIENT_MEMORY );
-        }
-
-    }
-    else
-    {
-        //TODO: reallocate
-    }
-
-    memcpy( workspace->H_EEM->start, Hptr->start, sizeof(unsigned int) * (system->N + 1) );
-    memcpy( workspace->H_EEM->j, Hptr->j, sizeof(unsigned int) * Hptr->start[system->N] );
-    memcpy( workspace->H_EEM->val, Hptr->val, sizeof(real) * Hptr->start[system->N] );
-    Hptr = workspace->H_EEM;
-
     switch ( control->cm_solver_pre_comp_type )
     {
     case DIAG_PC:
@@ -2130,9 +2202,7 @@ static void Setup_Preconditioner_EEM( const reax_system * const system,
 
         if ( workspace->L == NULL )
         {
-            if ( Allocate_Matrix( &(workspace->L_EEM), far_nbrs->n, fillin ) == FAILURE ||
-                    Allocate_Matrix( &(workspace->U_EEM), far_nbrs->n, fillin ) == FAILURE ||
-                    Allocate_Matrix( &(workspace->L), system->N_cm, fillin + system->N_cm ) == FAILURE ||
+            if ( Allocate_Matrix( &(workspace->L), system->N_cm, fillin + system->N_cm ) == FAILURE ||
                     Allocate_Matrix( &(workspace->U), system->N_cm, fillin + system->N_cm ) == FAILURE )
             {
                 fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
@@ -2150,9 +2220,7 @@ static void Setup_Preconditioner_EEM( const reax_system * const system,
         if ( workspace->L == NULL )
         {
             /* factors have sparsity pattern as H */
-            if ( Allocate_Matrix( &(workspace->L_EEM), system->N, Hptr->m ) == FAILURE ||
-                    Allocate_Matrix( &(workspace->U_EEM), system->N, Hptr->m ) == FAILURE ||
-                    Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
+            if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
                     Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
             {
                 fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
@@ -2175,9 +2243,7 @@ static void Setup_Preconditioner_EEM( const reax_system * const system,
         if ( workspace->L == NULL )
         {
             /* TODO: safest storage estimate is ILU(0) (same as lower triangular portion of H), could improve later */
-            if ( Allocate_Matrix( &(workspace->L_EEM), system->N, Hptr->m ) == FAILURE ||
-                    Allocate_Matrix( &(workspace->U_EEM), system->N, Hptr->m ) == FAILURE ||
-                    Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
+            if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
                     Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
             {
                 fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
@@ -2194,9 +2260,7 @@ static void Setup_Preconditioner_EEM( const reax_system * const system,
         if ( workspace->L == NULL )
         {
             /* factors have sparsity pattern as H */
-            if ( Allocate_Matrix( &(workspace->L_EEM), system->N, Hptr->m ) == FAILURE ||
-                    Allocate_Matrix( &(workspace->U_EEM), system->N, Hptr->m ) == FAILURE ||
-                    Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
+            if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
                     Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
             {
                 fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
@@ -2214,6 +2278,8 @@ static void Setup_Preconditioner_EEM( const reax_system * const system,
         exit( INVALID_INPUT );
         break;
     }
+
+    Hptr->val[Hptr->start[system->N + 1] - 1] = 0.0;
 }
 
 
@@ -2244,25 +2310,11 @@ static void Setup_Preconditioner_ACKS2( const reax_system * const system,
     {
         Sort_Matrix_Rows( workspace->H_sp );
     }
+    data->timing.cm_sort_mat_rows += Get_Timing_Info( time );
 
     Hptr->val[Hptr->start[system->N + 1] - 1] = 1.0;
     Hptr->val[Hptr->start[system->N_cm] - 1] = 1.0;
 
-    if ( control->cm_solver_pre_app_type == TRI_SOLVE_GC_PA )
-    {
-        if ( control->cm_domain_sparsify_enabled == TRUE )
-        {
-            Hptr = setup_graph_coloring( workspace->H_sp );
-        }
-        else
-        {
-            Hptr = setup_graph_coloring( workspace->H );
-        }
-
-        Sort_Matrix_Rows( Hptr );
-    }
-    data->timing.cm_sort_mat_rows += Get_Timing_Info( time );
-
 #if defined(DEBUG)
     fprintf( stderr, "H matrix sorted\n" );
 #endif
@@ -2397,13 +2449,19 @@ static void Calculate_Charges_QEq( const reax_system * const system,
     for ( i = 0; i < system->N_cm; ++i )
     {
         system->atoms[i].q = workspace->s[0][i] - u * workspace->t[0][i];
+
+#if defined(DEBUG_FOCUS)
+        printf("atom %4d: %f\n", i, system->atoms[i].q);
+        printf("  x[0]: %10.5f, x[1]: %10.5f, x[2]:  %10.5f\n",
+                system->atoms[i].x[0], system->atoms[i].x[1], system->atoms[i].x[2]);
+#endif
     }
 }
 
 
-/* Get atomic charge q for EEM method
+/* Get atomic charge q for EE method
  */
-static void Calculate_Charges_EEM( const reax_system * const system,
+static void Calculate_Charges_EE( const reax_system * const system,
         static_storage * const workspace )
 {
     int i;
@@ -2441,11 +2499,6 @@ static void QEq( reax_system * const system, control_params * const control,
 
     Extrapolate_Charges_QEq( system, control, data, workspace );
 
-//    if( data->step == 0 || data->step == 100 )
-//    {
-//      Print_Linear_System( system, control, workspace, data->step );
-//    }
-
     switch ( control->cm_solver_type )
     {
     case GMRES_S:
@@ -2507,17 +2560,20 @@ static void QEq( reax_system * const system, control_params * const control,
 
     Calculate_Charges_QEq( system, workspace );
 
-    //fprintf( stderr, "%d %.9f %.9f %.9f %.9f %.9f %.9f\n",
-    //   data->step,
-    //   workspace->s[0][0], workspace->t[0][0],
-    //   workspace->s[0][1], workspace->t[0][1],
-    //   workspace->s[0][2], workspace->t[0][2] );
-    // if( data->step == control->nsteps )
-    //Print_Charges( system, control, workspace, data->step );
+#if defined(DEBUG_FOCUS)
+    fprintf( stderr, "%d %.9f %.9f %.9f %.9f %.9f %.9f\n", data->step,
+       workspace->s[0][0], workspace->t[0][0],
+       workspace->s[0][1], workspace->t[0][1],
+       workspace->s[0][2], workspace->t[0][2] );
+    if( data->step == control->nsteps )
+    {
+        Print_Charges( system, control, workspace, data->step );
+    }
+#endif
 }
 
 
-/* Main driver method for EEM kernel
+/* Main driver method for EE kernel
  *
  * Rough outline:
  *  1) init / setup routines for preconditioning of linear solver
@@ -2526,7 +2582,7 @@ static void QEq( reax_system * const system, control_params * const control,
  *  4) perform 1 linear solve
  *  5) compute atomic charges based on output of (4)
  */
-static void EEM( reax_system * const system, control_params * const control,
+static void EE( reax_system * const system, control_params * const control,
         simulation_data * const data, static_storage * const workspace,
         const list * const far_nbrs, const output_controls * const out_control )
 {
@@ -2535,14 +2591,12 @@ static void EEM( reax_system * const system, control_params * const control,
     if ( control->cm_solver_pre_comp_refactor > 0 &&
             ((data->step - data->prev_steps) % control->cm_solver_pre_comp_refactor == 0) )
     {
-        Setup_Preconditioner_EEM( system, control, data, workspace, far_nbrs );
+        Setup_Preconditioner_EE( system, control, data, workspace, far_nbrs );
 
-        Compute_Preconditioner_EEM( system, control, data, workspace, far_nbrs );
+        Compute_Preconditioner_EE( system, control, data, workspace, far_nbrs );
     }
 
-//   Print_Linear_System( system, control, workspace, data->step );
-
-    Extrapolate_Charges_EEM( system, control, data, workspace );
+    Extrapolate_Charges_EE( system, control, data, workspace );
 
     switch ( control->cm_solver_type )
     {
@@ -2571,7 +2625,7 @@ static void EEM( reax_system * const system, control_params * const control,
         break;
 
     default:
-        fprintf( stderr, "Unrecognized QEq solver selection. Terminating...\n" );
+        fprintf( stderr, "Unrecognized EE solver selection. Terminating...\n" );
         exit( INVALID_INPUT );
         break;
     }
@@ -2582,7 +2636,7 @@ static void EEM( reax_system * const system, control_params * const control,
     fprintf( stderr, "linsolve-" );
 #endif
 
-    Calculate_Charges_EEM( system, workspace );
+    Calculate_Charges_EE( system, workspace );
 
     // if( data->step == control->nsteps )
     //Print_Charges( system, control, workspace, data->step );
@@ -2614,7 +2668,7 @@ static void ACKS2( reax_system * const system, control_params * const control,
 
 //   Print_Linear_System( system, control, workspace, data->step );
 
-    Extrapolate_Charges_EEM( system, control, data, workspace );
+    Extrapolate_Charges_EE( system, control, data, workspace );
 
     switch ( control->cm_solver_type )
     {
@@ -2654,7 +2708,7 @@ static void ACKS2( reax_system * const system, control_params * const control,
     fprintf( stderr, "linsolve-" );
 #endif
 
-    Calculate_Charges_EEM( system, workspace );
+    Calculate_Charges_EE( system, workspace );
 }
 
 
@@ -2663,16 +2717,18 @@ void Compute_Charges( reax_system * const system, control_params * const control
                       const list * const far_nbrs, const output_controls * const out_control )
 {
     int i;
-//    char fname[200];
-//    FILE * fp;
+#if defined(DEBUG_FOCUS)
+    char fname[200];
+    FILE * fp;
 
-//    if ( data->step == 10 || data->step == 50 || data->step == 100 )
-//    {
-//        sprintf( fname, "s_%d_%s.out", data->step, control->sim_name );
-//        fp = fopen( fname, "w" );
-//        Vector_Print( fp, NULL, workspace->s[0], system->N_cm );
-//        fclose( fp );
-//    }
+    if ( data->step == 0 || data->step == 10 || data->step == 50 || data->step == 100 )
+    {
+        sprintf( fname, "s_%d_%s.out", data->step, control->sim_name );
+        fp = fopen( fname, "w" );
+        Vector_Print( fp, NULL, workspace->s[0], system->N_cm );
+        fclose( fp );
+    }
+#endif
 
     switch ( control->charge_method )
     {
@@ -2680,8 +2736,8 @@ void Compute_Charges( reax_system * const system, control_params * const control
         QEq( system, control, data, workspace, far_nbrs, out_control );
         break;
 
-    case EEM_CM:
-        EEM( system, control, data, workspace, far_nbrs, out_control );
+    case EE_CM:
+        EE( system, control, data, workspace, far_nbrs, out_control );
         break;
 
     case ACKS2_CM:
@@ -2694,13 +2750,20 @@ void Compute_Charges( reax_system * const system, control_params * const control
         break;
     }
 
-//    if ( data->step == 0 )
-//    {
-//        sprintf( fname, "H_%d_%s.out", data->step, control->sim_name );
-//        Print_Sparse_Matrix2( workspace->H, fname );
-//
-//        sprintf( fname, "b_%d_%s.out", data->step, control->sim_name );
-//        fp = fopen( fname, "w" );
-//        Vector_Print( fp, NULL, workspace->b_s, system->N_cm );
-//    }
+#if defined(DEBUG_FOCUS)
+    if ( data->step == 0 || data->step == 10 || data->step == 50 || data->step == 100 )
+    {
+        if ( data->step == 0 )
+        {
+            sprintf( fname, "H_%d_%s.out", data->step, control->sim_name );
+            Print_Sparse_Matrix2( workspace->H, fname, NULL );
+            Print_Sparse_Matrix_Binary( workspace->H, fname );
+        }
+
+        sprintf( fname, "b_%d_%s.out", data->step, control->sim_name );
+        fp = fopen( fname, "w" );
+        Vector_Print( fp, NULL, workspace->b_s, system->N_cm );
+        fclose( fp );
+    }
+#endif
 }
diff --git a/sPuReMD/src/forces.c b/sPuReMD/src/forces.c
index e6c82684a4ce4a5dcbc6b21d772cdd3b98015bfb..2d5e9a49e87d64b1dd0e4fdd439c9bb51f177c33 100644
--- a/sPuReMD/src/forces.c
+++ b/sPuReMD/src/forces.c
@@ -310,7 +310,8 @@ static inline real Init_Charge_Matrix_Entry_Tab( reax_system *system,
         }
         break;
 
-    case EEM_CM:
+    case EE_CM:
+        //TODO
         switch ( pos )
         {
             case OFF_DIAGONAL:
@@ -388,7 +389,7 @@ static inline real Init_Charge_Matrix_Entry( reax_system *system,
         }
         break;
 
-    case EEM_CM:
+    case EE_CM:
         switch ( pos )
         {
             case OFF_DIAGONAL:
@@ -481,7 +482,7 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system *system,
         case QEQ_CM:
             break;
 
-        case EEM_CM:
+        case EE_CM:
             H->start[system->N_cm - 1] = *Htop;
             H_sp->start[system->N_cm - 1] = *H_sp_top;
 
@@ -1203,7 +1204,7 @@ void Init_Forces_Tab( reax_system *system, control_params *control,
              data->step, Htop, num_bonds, num_hbonds );
     //Print_Bonds( system, bonds, "sbonds.out" );
     //Print_Bond_List2( system, bonds, "sbonds.out" );
-    //Print_Sparse_Matrix2( H, "H.out" );
+    //Print_Sparse_Matrix2( H, "H.out", NULL );
 #endif
 }
 
diff --git a/sPuReMD/src/init_md.c b/sPuReMD/src/init_md.c
index 4b93ec83c4d6dcfbb8c32ba2c7736889ba6fc635..4228f8bd4967cca1822e5d07fd0d5055503d5b20 100644
--- a/sPuReMD/src/init_md.c
+++ b/sPuReMD/src/init_md.c
@@ -298,7 +298,7 @@ void Init_Workspace( reax_system *system, control_params *control,
         case QEQ_CM:
             system->N_cm = system->N;
             break;
-        case EEM_CM:
+        case EE_CM:
             system->N_cm = system->N + 1;
             break;
         case ACKS2_CM:
@@ -314,15 +314,6 @@ void Init_Workspace( reax_system *system, control_params *control,
     workspace->H_sp = NULL;
     workspace->L = NULL;
     workspace->U = NULL;
-    workspace->H_EEM = NULL;
-    workspace->L_EEM = NULL;
-    workspace->U_EEM = NULL;
-    workspace->H_ACKS2_1 = NULL;
-    workspace->H_ACKS2_2 = NULL;
-    workspace->L_ACKS2_1 = NULL;
-    workspace->L_ACKS2_2 = NULL;
-    workspace->U_ACKS2_1 = NULL;
-    workspace->U_ACKS2_2 = NULL;
     workspace->Hdia_inv = NULL;
     if ( control->cm_solver_pre_comp_type == ICHOLT_PC ||
             control->cm_solver_pre_comp_type == ILUT_PAR_PC )
@@ -361,7 +352,7 @@ void Init_Workspace( reax_system *system, control_params *control,
             }
             break;
 
-        case EEM_CM:
+        case EE_CM:
             for ( i = 0; i < system->N; ++i )
             {
                 workspace->b_s[i] = -system->reaxprm.sbp[ system->atoms[i].type ].chi;
@@ -535,7 +526,7 @@ void Init_Lists( reax_system *system, control_params *control,
         case QEQ_CM:
             max_nnz = Htop;
             break;
-        case EEM_CM:
+        case EE_CM:
             max_nnz = Htop + system->N_cm;
             break;
         case ACKS2_CM:
diff --git a/sPuReMD/src/mytypes.h b/sPuReMD/src/mytypes.h
index 7487086ec260f15f8731834ca444c02638365d17..b1090a2527aa879c00b1caa05c49839746bf87cf 100644
--- a/sPuReMD/src/mytypes.h
+++ b/sPuReMD/src/mytypes.h
@@ -193,7 +193,7 @@ enum geo_formats
 
 enum charge_method
 {
-    QEQ_CM = 0, EEM_CM = 1, ACKS2_CM = 2,
+    QEQ_CM = 0, EE_CM = 1, ACKS2_CM = 2,
 };
 
 enum solver
@@ -781,10 +781,6 @@ typedef struct
 
     /* charge method storage */
     sparse_matrix *H, *H_sp, *L, *U;
-    /* EEM-specific */
-    sparse_matrix *H_EEM, *L_EEM, *U_EEM;
-    /* ACKS2-specific */
-    sparse_matrix *H_ACKS2_1, *H_ACKS2_2, *L_ACKS2_1, *L_ACKS2_2, *U_ACKS2_1, *U_ACKS2_2;
     real *droptol;
     real *w;
     real *Hdia_inv;
diff --git a/sPuReMD/src/print_utils.c b/sPuReMD/src/print_utils.c
index 8a5952760cd10372aa0ab548e0ed4886e85abfab..06b7d63dfbf0420a1b062ebc39a6c5bfee0f9234 100644
--- a/sPuReMD/src/print_utils.c
+++ b/sPuReMD/src/print_utils.c
@@ -827,10 +827,19 @@ void Print_Sparse_Matrix( sparse_matrix *A )
 }
 
 
-void Print_Sparse_Matrix2( sparse_matrix *A, char *fname )
+void Print_Sparse_Matrix2( sparse_matrix *A, char *fname, char *mode )
 {
     int i, j;
-    FILE *f = fopen( fname, "w" );
+    FILE *f;
+   
+    if ( mode == NULL )
+    {
+        f = fopen( fname, "w" );
+    }
+    else
+    {
+        f = fopen( fname, mode );
+    }
 
     for ( i = 0; i < A->n; ++i )
     {
@@ -846,6 +855,43 @@ void Print_Sparse_Matrix2( sparse_matrix *A, char *fname )
 }
 
 
+/* Note: watch out for portability issues with endianness
+ * due to serialization of numeric types (integer, IEEE 754) */
+void Print_Sparse_Matrix_Binary( sparse_matrix *A, char *fname )
+{
+    int i, j, temp;
+    FILE *f;
+   
+    f = fopen( fname, "wb" );
+
+    /* header: # rows, # nonzeros */
+    fwrite( &(A->n), sizeof(unsigned int), 1, f );
+    fwrite( &(A->start[A->n]), sizeof(unsigned int), 1, f );
+
+    /* row pointers */
+    for ( i = 0; i <= A->n; ++i )
+    {
+        //Convert 0-based to 1-based (for Matlab)
+        temp = A->start[i] + 1;
+        fwrite( &temp, sizeof(unsigned int), 1, f );
+    }
+
+    /* column indices and non-zeros */
+    for ( i = 0; i <= A->n; ++i )
+    {
+        for ( j = A->start[i]; j < A->start[i + 1]; ++j )
+        {
+            //Convert 0-based to 1-based (for Matlab)
+            temp = A->j[j] + 1;
+            fwrite( &temp, sizeof(unsigned int), 1, f );
+            fwrite( &(A->val[j]), sizeof(real), 1, f );
+        }
+    }
+
+    fclose(f);
+}
+
+
 void Print_Bonds( reax_system *system, list *bonds, char *fname )
 {
     int i, pj;
diff --git a/sPuReMD/src/print_utils.h b/sPuReMD/src/print_utils.h
index 8b3b363571ca5b7c50195b4c91bca8951c1caa7b..8c15ebbe81487eabf2a5db13627eef60ee8432d3 100644
--- a/sPuReMD/src/print_utils.h
+++ b/sPuReMD/src/print_utils.h
@@ -55,7 +55,9 @@ void Print_Soln( static_storage*, real*, real*, real*, int );
 
 void Print_Sparse_Matrix( sparse_matrix* );
 
-void Print_Sparse_Matrix2( sparse_matrix*, char* );
+void Print_Sparse_Matrix2( sparse_matrix*, char*, char* );
+
+void Print_Sparse_Matrix_Binary( sparse_matrix*, char* );
 
 void Print_Bonds( reax_system*, list*, char* );
 void Print_Bond_List2( reax_system*, list*, char* );
diff --git a/sPuReMD/src/vector.c b/sPuReMD/src/vector.c
index 483ba806d45432f9e2e0a0029d1db65058229cad..5cabd03b7f055ad02607eea052f5988c6ed2eaa3 100644
--- a/sPuReMD/src/vector.c
+++ b/sPuReMD/src/vector.c
@@ -114,7 +114,7 @@ inline void Vector_Add( real * const dest, const real c, const real * const v, c
 
 
 void Vector_Print( FILE * const fout, const char * const vname, const real * const v,
-                   const unsigned int k )
+        const unsigned int k )
 {
     unsigned int i;