MOLINFO WHOLE: check that the reference is not split by the periodic boundary (#1492)

CHANGES/v2.11.md

@@ -5,6 +5,7 @@ This page contains changes that will end up in 2.11
 ### Changes relevant for users:
 
   - When using [MOLINFO](MOLINFO.md) with the `WHOLE` flag, PBCs in the following actions will be reconstructed using a minimum spanning tree based on the coordinates stored in the MOLFILE reference pdb.
+  - When the reference pdb of [MOLINFO](MOLINFO.md) with the `WHOLE` flag has a CRYST1 record with space group P 1, PLUMED now checks along each chain that the reference is whole in that box, and writes a warning to the log naming the atoms that are across the periodic boundary from the atom before them. Before, a split reference made [WHOLEMOLECULES](WHOLEMOLECULES.md) split whole molecules without any warning.
   - When using the benchmark now it is possible to scale the atom distances in the synthetic atom distribution (the base atomic distance is 1)
   - [SAXS](SAXS.md) now accepts the `MARTINI_FF` keyword to select the Martini form-factor generation: `MARTINI_FF=2` (default) keeps the existing Martini 2.2 protein and nucleic-acid form factors, while `MARTINI_FF=3` selects Martini 3 protein form factors.
   - Added environment variable `PLUMED_LOG_FILE` to control the log file pathname when runnning with GROMACS (the log file was hardcoded as `PLUMED.OUT`)

regtest/basic/rt-molinfo-whole-check/Makefile

@@ -0,0 +1 @@
+include ../../scripts/test.make

regtest/basic/rt-molinfo-whole-check/config

@@ -0,0 +1 @@
+type=make

regtest/basic/rt-molinfo-whole-check/main.cpp

@@ -0,0 +1,62 @@
+#include "plumed/wrapper/Plumed.h"
+#include <fstream>
+#include <sstream>
+#include <string>
+#include <vector>
+
+// MOLINFO WHOLE takes the reference structure as whole. Here each reference is read with
+// MOLINFO ... WHOLE and used by WHOLEMOLECULES, and we report whether PLUMED accepts it or
+// writes a warning that a chain looks split by the periodic boundary of its CRYST1 record.
+// torn.pdb is whole.pdb with 15 atoms wrapped across a face of its (rhombic dodecahedron) cell.
+// Its copies without a usable cell (no CRYST1, the 1 A placeholder, a zero cell, a crystallographic
+// space group) cannot be checked and are accepted.
+// dna_torn.pdb is dna_whole.pdb (chain A of 1BNA) wrapped nucleotide by nucleotide: each nucleotide is
+// whole, and the chain is split between them.
+int main(int, char**) {
+  std::ofstream out("output");
+  unsigned natoms=304;
+  std::vector<double> positions(3*natoms,0.0);
+  const std::vector<std::pair<std::string,std::string>> cases{
+    {"whole.pdb","1-304"},
+    {"torn.pdb","1-304"},
+    {"torn_nocryst.pdb","1-304"},
+    {"torn_placeholder.pdb","1-304"},
+    {"torn_zero.pdb","1-304"},
+    {"torn_crystal.pdb","1-304"},
+    {"dna_whole.pdb","1-243"},
+    {"dna_torn.pdb","1-243"}
+  };
+  for(unsigned i=0; i<cases.size(); ++i) {
+    const auto & c(cases[i]);
+    const std::string logfile="test"+std::to_string(i)+".log";
+    out<<c.first<<" ENTITY0="<<c.second<<" : ";
+    try {
+      PLMD::Plumed plumed;
+      plumed.cmd("setNatoms",&natoms);
+      plumed.cmd("setLogFile",logfile.c_str());
+      plumed.cmd("init");
+      plumed.cmd("readInputLine",("MOLINFO STRUCTURE="+c.first+" WHOLE").c_str());
+      plumed.cmd("readInputLine",("WHOLEMOLECULES ENTITY0="+c.second).c_str());
+    } catch(PLMD::Plumed::ExceptionError &e) {
+      out<<"unexpected error: "<<e.what()<<"\n";
+      continue;
+    }
+    // the log is complete once plumed is destroyed
+    std::ifstream in(logfile);
+    std::stringstream ss;
+    ss<<in.rdbuf();
+    const std::string m(ss.str());
+    const auto a=m.find("may not be whole: ");
+    if(a==std::string::npos) {
+      out<<"accepted\n";
+      continue;
+    }
+    const auto e=m.find(". WHOLEMOLECULES",a);
+    if(e==std::string::npos) {
+      out<<"unexpected warning: "<<m.substr(a)<<"\n";
+    } else {
+      out<<"warning, "<<m.substr(a+18,e-a-18)<<"\n";
+    }
+  }
+  return 0;
+}

regtest/basic/rt-molinfo-whole-check/output.reference

@@ -0,0 +1,8 @@
+whole.pdb ENTITY0=1-304 : accepted
+torn.pdb ENTITY0=1-304 : warning, in the box of its CRYST1 record, 14 atoms are across the periodic boundary from the atom before them in their chain: 200 (SER13) 202 (SER13) 205 (SER13) 206 (SER13) 207 (SER13) 208 (SER13) 211 (SER14) 213 (SER14) 216 (SER14) 217 (SER14) and 4 more
+torn_nocryst.pdb ENTITY0=1-304 : accepted
+torn_placeholder.pdb ENTITY0=1-304 : accepted
+torn_zero.pdb ENTITY0=1-304 : accepted
+torn_crystal.pdb ENTITY0=1-304 : accepted
+dna_whole.pdb ENTITY0=1-243 : accepted
+dna_torn.pdb ENTITY0=1-243 : warning, in the box of its CRYST1 record, 1 atoms are across the periodic boundary from the atom before them in their chain: 122 (DT7)