Do not set dummy unit cell vectors along non-periodic dimensions 34/4334/2
authorGiacomo Fiorin <giacomo.fiorin@gmail.com>
Sat, 7 Jul 2018 15:18:47 +0000 (11:18 -0400)
committerGiacomo Fiorin <giacomo.fiorin@gmail.com>
Sat, 14 Jul 2018 21:05:44 +0000 (17:05 -0400)
If no periodicity is defined along a certain dimension, NAMD currently sets a
"dummy" PBC unit cell vector (e.g. (1,0,0) for the first dimension).

Due to this choice, analysis tools cannot properly detect the boundary
conditions of the simulation from DCD files, which do not support
integer/boolean periodicity flags.  The unit cells stored in XST and XSC files
are correct (non-periodic dimensions are omitted), but these files are very
rarely used in analysis.

The above behavior was likely introduced to avoid divisions by zero when
computing reciprocal lattice vectors, but it can potentially cause undefined
behavior in various components of NAMD unless full PBCs are enabled.

This change partially addresses the problem by setting non-periodic dimensions
to zero when writing a DCD file.

The desirable solution is to eliminate all uses of the "dummy" unit cell
vectors altogether.

Change-Id: I9983b070b2a5f161d5d9859a72690b955d574077

src/Output.C

index f07276c..7da2987 100644 (file)
@@ -129,26 +129,30 @@ void NAMD_close(int fd, const char *fname) {
 
 
 static void lattice_to_unitcell(const Lattice *lattice, double *unitcell) {
-   if (lattice && lattice->a_p() && lattice->b_p() && lattice->c_p()) {
-      const Vector &a=lattice->a();
-      const Vector &b=lattice->b();
-      const Vector &c=lattice->c();
-      unitcell[0] = a.length();
-      unitcell[2] = b.length();
-      unitcell[5] = c.length();
-      double cosAB = (a*b)/(unitcell[0]*unitcell[2]);
-      double cosAC = (a*c)/(unitcell[0]*unitcell[5]);
-      double cosBC = (b*c)/(unitcell[2]*unitcell[5]);
-      if (cosAB > 1.0) cosAB = 1.0; else if (cosAB < -1.0) cosAB = -1.0;
-      if (cosAC > 1.0) cosAC = 1.0; else if (cosAC < -1.0) cosAC = -1.0;
-      if (cosBC > 1.0) cosBC = 1.0; else if (cosBC < -1.0) cosBC = -1.0;
-      unitcell[1] = cosAB;
-      unitcell[3] = cosAC;
-      unitcell[4] = cosBC;
-   } else {
-      unitcell[0] = unitcell[2] = unitcell[5] = 1.0;
-      unitcell[1] = unitcell[3] = unitcell[4] = 0.0;
-   }
+
+  unitcell[0] = unitcell[2] = unitcell[5] = 0.0;
+  unitcell[1] = unitcell[3] = unitcell[4] = 0.0;
+
+  if (lattice) {
+    const Vector &a=lattice->a();
+    const Vector &b=lattice->b();
+    const Vector &c=lattice->c();
+    unitcell[0] = (lattice->a_p()) ? a.length() : 0.0;
+    unitcell[2] = (lattice->b_p()) ? b.length() : 0.0;
+    unitcell[5] = (lattice->c_p()) ? c.length() : 0.0;
+    double cosAB = (lattice->a_p() && lattice->b_p() ) ?
+      (a*b)/(unitcell[0]*unitcell[2]) : 0.0;
+    double cosAC = (lattice->a_p() && lattice->c_p() ) ?
+      (a*c)/(unitcell[0]*unitcell[5]) : 0.0;
+    double cosBC = (lattice->b_p() && lattice->c_p() ) ?
+      (b*c)/(unitcell[2]*unitcell[5]) : 0.0;
+    if (cosAB > 1.0) cosAB = 1.0; else if (cosAB < -1.0) cosAB = -1.0;
+    if (cosAC > 1.0) cosAC = 1.0; else if (cosAC < -1.0) cosAC = -1.0;
+    if (cosBC > 1.0) cosBC = 1.0; else if (cosBC < -1.0) cosBC = -1.0;
+    unitcell[1] = cosAB;
+    unitcell[3] = cosAC;
+    unitcell[4] = cosBC;
+  }
 }