Add DSSP v4 support to do_dssp

[alexxy/gromacs.git] / src / gromacs / fileio / pdbio.cpp

diff --git a/src/gromacs/fileio/pdbio.cpp b/src/gromacs/fileio/pdbio.cpp
index 2e5d1ec39f89b2eae176300556a9ff94ae49cf17..c3c3995039774070b923690aa6330de820617a77 100644 (file)
--- a/src/gromacs/fileio/pdbio.cpp
+++ b/src/gromacs/fileio/pdbio.cpp
@@ -49,11 +49,11 @@
 
 #include "gromacs/fileio/gmxfio.h"
 #include "gromacs/math/units.h"
+#include "gromacs/math/utilities.h"
 #include "gromacs/math/vec.h"
 #include "gromacs/pbcutil/pbc.h"
 #include "gromacs/topology/atomprop.h"
 #include "gromacs/topology/ifunc.h"
-#include "gromacs/topology/residuetypes.h"
 #include "gromacs/topology/symtab.h"
 #include "gromacs/topology/topology.h"
 #include "gromacs/utility/coolstuff.h"
@@ -86,8 +86,6 @@ const char* enumValueToString(PdbRecordType enumValue)
     return pdbRecordTypeName[enumValue];
 }
 
-#define REMARK_SIM_BOX "REMARK    THIS IS A SIMULATION BOX"
-
 void gmx_write_pdb_box(FILE* out, PbcType pbcType, const matrix box)
 {
     real alpha, beta, gamma;
@@ -104,7 +102,7 @@ void gmx_write_pdb_box(FILE* out, PbcType pbcType, const matrix box)
 
     if (norm2(box[YY]) * norm2(box[ZZ]) != 0)
     {
-        alpha = RAD2DEG * gmx_angle(box[YY], box[ZZ]);
+        alpha = gmx::c_rad2Deg * gmx_angle(box[YY], box[ZZ]);
     }
     else
     {
@@ -112,7 +110,7 @@ void gmx_write_pdb_box(FILE* out, PbcType pbcType, const matrix box)
     }
     if (norm2(box[XX]) * norm2(box[ZZ]) != 0)
     {
-        beta = RAD2DEG * gmx_angle(box[XX], box[ZZ]);
+        beta = gmx::c_rad2Deg * gmx_angle(box[XX], box[ZZ]);
     }
     else
     {
@@ -120,7 +118,7 @@ void gmx_write_pdb_box(FILE* out, PbcType pbcType, const matrix box)
     }
     if (norm2(box[XX]) * norm2(box[YY]) != 0)
     {
-        gamma = RAD2DEG * gmx_angle(box[XX], box[YY]);
+        gamma = gmx::c_rad2Deg * gmx_angle(box[XX], box[YY]);
     }
     else
     {
@@ -206,7 +204,7 @@ static void read_cryst1(char* line, PbcType* pbcType, matrix box)
         {
             if (alpha != 90.0)
             {
-                cosa = std::cos(alpha * DEG2RAD);
+                cosa = std::cos(alpha * gmx::c_deg2Rad);
             }
             else
             {
@@ -214,7 +212,7 @@ static void read_cryst1(char* line, PbcType* pbcType, matrix box)
             }
             if (beta != 90.0)
             {
-                cosb = std::cos(beta * DEG2RAD);
+                cosb = std::cos(beta * gmx::c_deg2Rad);
             }
             else
             {
@@ -222,8 +220,8 @@ static void read_cryst1(char* line, PbcType* pbcType, matrix box)
             }
             if (gamma != 90.0)
             {
-                cosg = std::cos(gamma * DEG2RAD);
-                sing = std::sin(gamma * DEG2RAD);
+                cosg = std::cos(gamma * gmx::c_deg2Rad);
+                sing = std::sin(gamma * gmx::c_deg2Rad);
             }
             else
             {
@@ -298,7 +296,8 @@ void write_pdbfile_indexed(FILE*          out,
                            int            nindex,
                            const int      index[],
                            gmx_conect     conect,
-                           bool           usePqrFormat)
+                           bool           usePqrFormat,
+                           bool           standardCompliantMode)
 {
     gmx_conect_t* gc = static_cast<gmx_conect_t*>(conect);
     PdbRecordType type;
@@ -306,8 +305,27 @@ void write_pdbfile_indexed(FILE*          out,
     real          occup, bfac;
     gmx_bool      bOccup;
 
+    if (standardCompliantMode)
+    {
+        fprintf(out, "HEADER    GROMACS SIMULATION BOX                  01-JAN-00   0000\n");
+        fprintf(out,
+                "TITLE     Gromacs simulation box\n"
+                "COMPND    MOL_ID:  1;                                                           \n"
+                "COMPND   2 MOLECULE:  GROMACS SIMULATION BOX;                                   \n"
+                "COMPND   3 CHAIN: A;  \n"
+                "SOURCE    MOL_ID: 1;\n"
+                "SOURCE   2 SYNTHETIC\n"
+                "KEYWDS    GROMACS\n"
+                "EXPDTA    PURE PRODUCT OF COMPUTER SIMULATION\n"
+                "AUTHOR    GROMACS\n"
+                "REVDAT   1   01-JAN-00 0000    0\n");
+    }
+    else
+    {
+        fprintf(out, "TITLE     %s\n", (title && title[0]) ? title : gmx::bromacs().c_str());
+    }
+
 
-    fprintf(out, "TITLE     %s\n", (title && title[0]) ? title : gmx::bromacs().c_str());
     if (box && ((norm2(box[XX]) != 0.0F) || (norm2(box[YY]) != 0.0F) || (norm2(box[ZZ]) != 0.0F)))
     {
         gmx_write_pdb_box(out, pbcType, box);
@@ -331,7 +349,11 @@ void write_pdbfile_indexed(FILE*          out,
 
     fprintf(out, "MODEL %8d\n", model_nr > 0 ? model_nr : 1);
 
-    ResidueType rt;
+    // Collect last printed values for TER record
+    int         lastAtomNumber = 0;
+    std::string lastResName;
+    int         lastResNr = 0;
+
     for (int ii = 0; ii < nindex; ii++)
     {
         int         i      = index[ii];
@@ -394,6 +416,10 @@ void write_pdbfile_indexed(FILE*          out,
                                      bfac,
                                      atoms->atom[i].elem);
 
+            lastAtomNumber = i + 1;
+            lastResName    = resnm;
+            lastResNr      = resnr;
+
             if (atoms->pdbinfo && atoms->pdbinfo[i].bAnisotropic)
             {
                 fprintf(out,
@@ -429,7 +455,15 @@ void write_pdbfile_indexed(FILE*          out,
         }
     }
 
-    fprintf(out, "TER\n");
+    if (standardCompliantMode)
+    {
+        fprintf(out, "TER   %5d      %4.4s%c%4d\n", lastAtomNumber, lastResName.c_str(), chainid, lastResNr);
+    }
+    else
+    {
+        fprintf(out, "TER\n");
+    }
+
     fprintf(out, "ENDMDL\n");
 
     if (nullptr != gc)
@@ -574,18 +608,20 @@ void get_pdb_atomnumber(const t_atoms* atoms, AtomProperties* aps)
                 atomNumberSet = true;
             }
         }
-        std::string buf;
+        static constexpr size_t sc_maxElementNameLength = 3;
+        static_assert(sizeof(atoms->atom[i].elem) >= sc_maxElementNameLength + 1);
+        std::string element;
         if (atomNumberSet)
         {
             atoms->atom[i].atomnumber = atomnumber;
-            buf                       = aps->elementFromAtomNumber(atomnumber);
+            element                   = aps->elementFromAtomNumber(atomnumber);
             if (debug)
             {
                 fprintf(debug, "Atomnumber for atom '%s' is %d\n", anm, atomnumber);
             }
         }
-        buf.resize(3);
-        std::strncpy(atoms->atom[i].elem, buf.c_str(), 4);
+        element.resize(sc_maxElementNameLength);
+        std::strcpy(atoms->atom[i].elem, element.c_str());
     }
 }
 
@@ -873,13 +909,14 @@ int read_pdbfile(FILE*      in,
     title[0] = '\0';
     natom    = 0;
     chainnum = 0;
-    static const gmx::StringToEnumValueConverter<PdbRecordType, enumValueToString, gmx::StringCompareType::Exact, gmx::StripStrings::Yes> pdbRecordTypeIdentifier;
+    static const gmx::StringToEnumValueConverter<PdbRecordType, enumValueToString, gmx::StringCompareType::Exact, gmx::StripStrings::Yes>
+            sc_pdbRecordTypeIdentifier;
     while (!bStop && (fgets2(line, STRLEN, in) != nullptr))
     {
         // Convert line to a null-terminated string, then take a substring of at most 6 chars
         std::string                  lineAsString(line);
         std::optional<PdbRecordType> line_type =
-                pdbRecordTypeIdentifier.valueFrom(lineAsString.substr(0, 6));
+                sc_pdbRecordTypeIdentifier.valueFrom(lineAsString.substr(0, 6));
         if (!line_type)
         {
             // Skip this line because it does not start with a