Import/export of dcd format is implemented.

git-svn-id: svn+ssh://svn.sourceforge.jp/svnroot/molby/trunk@146 a2be9bc6-48de-4e38-9406-05402d4bc13c

diff --git a/MolLib/Dcd.c b/MolLib/Dcd.c
index a1fc3dc..b6a5007 100644 (file)
--- a/MolLib/Dcd.c
+++ b/MolLib/Dcd.c
@@ -25,6 +25,14 @@ s_Swap4(char *cp)
        cp[0] = w[3]; cp[1] = w[2]; cp[2] = w[1]; cp[3] = w[0];
 }
 
+static void
+s_Swap8(char *cp)
+{
+       char w[8];
+       w[0] = cp[0]; w[1] = cp[1]; w[2] = cp[2]; w[3] = cp[3]; w[4] = cp[4]; w[5] = cp[5]; w[6] = cp[6]; w[7] = cp[7];
+       cp[0] = w[7]; cp[1] = w[6]; cp[2] = w[5]; cp[3] = w[4]; cp[4] = w[3]; cp[5] = w[2]; cp[6] = w[1]; cp[7] = w[0];
+}
+
 #define s_SwapInt32(ip)  s_Swap4((char *)(ip))
 #define s_SwapSFloat32(fp) s_Swap4((char *)(fp))
 
@@ -48,6 +56,16 @@ s_Write4(int fd, const char *cp, int swap)
 }
 
 static int
+s_Write8(int fd, const char *cp, int swap)
+{
+       if (swap) {
+               char w[8];
+               w[7] = cp[0]; w[6] = cp[1]; w[5] = cp[2]; w[4] = cp[3]; w[3] = cp[4]; w[2] = cp[5]; w[1] = cp[6]; w[0] = cp[7];
+               return write(fd, w, 8);
+       } else return write(fd, cp, 8);
+}
+
+static int
 s_WriteInt32(DcdRecord *dcd, Int32 i)
 {
        return (s_Write4(dcd->fd, (const char *)(&i), dcd->reverse_endian) == 4);
@@ -60,6 +78,12 @@ s_WriteSFloat32(DcdRecord *dcd, SFloat32 f)
 }
 
 static int
+s_WriteSFloat64(DcdRecord *dcd, SFloat64 d)
+{
+       return (s_Write8(dcd->fd, (const char *)(&d), dcd->reverse_endian) == 8);
+}
+
+static int
 s_WriteZeros(int fd, int count)
 {
        char buf[128];
@@ -78,7 +102,8 @@ DcdOpen(const char *name, DcdRecord *dr)
 {
     Int32 nn, nnn;
     char buf[5];
-       
+       char delta_buf[8];
+
        if (dr == NULL)
                return -1;  /*  Internal error  */
        memset(dr, 0, sizeof(DcdRecord));
@@ -86,8 +111,9 @@ DcdOpen(const char *name, DcdRecord *dr)
        if (dr->fd < 0)
                return -2;  /*  Cannot open file  */
        
-    /*  Section 1: 'CORD', NFILE, NPRIV, NSAVC, NSTEP, 5*ZERO, DELTA,  */
-    /* WITH_UNITCELL, 4*ZERO, 24  */
+       /*  Section 1: 'CORD', NFILE, NPRIV, NSAVC, NSTEP, 4*ZERO, NFIX,  */
+       /* DELTA (4 bytes for charmm format, 8 bytes for x-plor format),  */
+       /* NEXTRA (charmm only), N4DIM, 6*ZERO, NCHARM, 24 */
        if (!s_ReadInt32(dr, &nn))
                return 1;   /*  Bad format: premature EOF  */
        if (nn != 84) {
@@ -102,16 +128,28 @@ DcdOpen(const char *name, DcdRecord *dr)
        s_ReadInt32(dr, &(dr->nstart));
        s_ReadInt32(dr, &(dr->ninterval));
        s_ReadInt32(dr, &(dr->nend));
-    lseek(dr->fd, 20, SEEK_CUR);  /*  Skip 5 zeros  */
-    s_ReadSFloat32(dr, &(dr->delta));
-       s_ReadInt32(dr, &(dr->with_unitcell));
-       lseek(dr->fd, 32, SEEK_CUR);  /*  Skip 8 zeros  */
-       s_ReadInt32(dr, &nn);  /*  This should be 24  */
-       if (nn != 24)
-               return 4;   /*  Bad format: bad end of first section  */
+    lseek(dr->fd, 16, SEEK_CUR);  /*  Skip 4 zeros  */
+       s_ReadInt32(dr, &(dr->nfix));
+       read(dr->fd, delta_buf, 8);   /*  If charmm format, then this is { float; int32; } otherwise this is a double */
+       s_ReadInt32(dr, &(dr->n4dim));
+       lseek(dr->fd, 28, SEEK_CUR);  /*  Skip 7 zeros  */
+       s_ReadInt32(dr, &(dr->ncharmver));
        s_ReadInt32(dr, &nn);  /*  This should be 84  */
        if (nn != 84)
                return 4;
+       if (dr->ncharmver == 0) {
+               if (dr->reverse_endian)
+                       s_Swap8(delta_buf);
+               dr->delta = *((SFloat64 *)delta_buf);
+               dr->nextra = 0;
+       } else {
+               if (dr->reverse_endian) {
+                       s_Swap4(delta_buf);
+                       s_Swap4(delta_buf + 4);
+               }
+               dr->delta = *((SFloat32 *)delta_buf);
+               dr->nextra = *((Int32 *)(delta_buf + 4));
+       }
        
     /*  Section 2: Title lines  */
     if (!s_ReadInt32(dr, &nn) || lseek(dr->fd, nn, SEEK_CUR) < 0 || !s_ReadInt32(dr, &nnn) || nn != nnn)
@@ -140,19 +178,21 @@ DcdCreate(const char *name, DcdRecord *dr)
                return -2;  /*  Cannot create file  */
        memset(buf, ' ', 80);
        
-    /*  Section 1: 'CORD', NFILE, NPRIV, NSAVC, NSTEP, 5*ZERO, DELTA,  */
-    /* WITH_UNITCELL, 8*ZERO, 24 (92 bytes) */
+       /*  Section 1: 'CORD', NFILE, NPRIV, NSAVC, NSTEP, 4*ZERO, NFIX,  */
+       /* DELTA (4 bytes for charmm format, 8 bytes for x-plor format),  */
+       /* NEXTRA (charmm only), N4DIM, 6*ZERO, NCHARM, 24 */
        if (!s_WriteInt32(dr, 84) ||
                write(dr->fd, "CORD", 4) != 4 ||
                !s_WriteInt32(dr, dr->nframes) ||
                !s_WriteInt32(dr, dr->nstart) ||
                !s_WriteInt32(dr, dr->ninterval) || 
                !s_WriteInt32(dr, dr->nend) ||
-               !s_WriteZeros(dr->fd, 20) ||
+               !s_WriteZeros(dr->fd, 16) ||
+               !s_WriteInt32(dr, dr->nfix) ||
                !s_WriteSFloat32(dr, dr->delta) ||
-               !s_WriteInt32(dr, dr->with_unitcell) ||
+               !s_WriteInt32(dr, dr->nextra) ||
                !s_WriteZeros(dr->fd, 32) ||
-               !s_WriteInt32(dr, 24) ||
+               !s_WriteInt32(dr, dr->ncharmver) ||
                !s_WriteInt32(dr, 84))
                return 1;   /*  Cannot write  */
        
@@ -189,20 +229,34 @@ int
 DcdReadFrame(DcdRecord *dr, int index, SFloat32 *xp, SFloat32 *yp, SFloat32 *zp, SFloat32 *cellp)
 {
     Int32 nn, nnn, i;
-    off_t block_size = 24 + (off_t)(dr->natoms * 12) + (dr->with_unitcell ? 56 : 0);
+    off_t block_size = (index > 0 ? dr->block_size : 24 + (off_t)(dr->natoms * 12) + (dr->nextra ? 56 : 0));
     lseek(dr->fd, dr->header_size + block_size * index, SEEK_SET);
-    if (dr->with_unitcell) {
-        if (!s_ReadInt32(dr, &nn) || nn != 48)
-            goto error;
-        if (cellp != NULL) {
-                       read(dr->fd, cellp, 48);
-                       for (i = 0; i < 6; i++)
-                               s_SwapSFloat32(cellp + i);
+       block_size = 0;
+       if (!s_ReadInt32(dr, &nn))
+               goto error;
+       if (nn == 48) {
+               if (dr->nextra) {
+                       SFloat64 mycell[6];
+                       if (nn != 48)
+                               goto error;
+                       read(dr->fd, mycell, 48);
+                       if (cellp != NULL) {
+                               for (i = 0; i < 6; i++) {
+                                       if (dr->reverse_endian)
+                                               s_Swap8((char *)(mycell + i));
+                                       cellp[i] = mycell[i];  /*  double -> float  */
+                               }
+                       }
+                       if (!s_ReadInt32(dr, &nn) || nn != 48)
+                               goto error;
+               } else {
+                       lseek(dr->fd, 52, SEEK_CUR);
                }
-        if (!s_ReadInt32(dr, &nn) || nn != 48)
-            goto error;
-    }
-    if (!s_ReadInt32(dr, &nn) || nn != dr->natoms * 4 ||
+               if (!s_ReadInt32(dr, &nn))
+                       goto error;
+               block_size = 56;
+       }
+    if (nn != dr->natoms * 4 ||
                read(dr->fd, xp, nn) < nn ||
                !s_ReadInt32(dr, &nnn) || nn != nnn)
         goto error;
@@ -226,6 +280,9 @@ DcdReadFrame(DcdRecord *dr, int index, SFloat32 *xp, SFloat32 *yp, SFloat32 *zp,
                for (i = 0; i < dr->natoms; i++)
                        s_SwapSFloat32(zp + i);
        }
+       block_size += 24 + (off_t)(dr->natoms * 12);
+       if (index == 0)
+               dr->block_size = block_size;
     return 0;
 error:
     return 1;
@@ -238,11 +295,11 @@ int
 DcdWriteFrame(DcdRecord *dr, int index, const SFloat32 *xp, const SFloat32 *yp, const SFloat32 *zp, const SFloat32 *cellp)
 {
        Int32 i;
-    if (dr->with_unitcell) {
+    if (dr->nextra) {
         if (!s_WriteInt32(dr, 48))
             goto error;
                for (i = 0; i < 6; i++)
-                       s_WriteSFloat32(dr, (cellp != NULL ? cellp[i] : 0.0));
+                       s_WriteSFloat64(dr, (cellp != NULL ? cellp[i] : 0.0));
                if (!s_WriteInt32(dr, 48))
             goto error;
     }

diff --git a/MolLib/Dcd.h b/MolLib/Dcd.h
index 495a9ce..2a3b474 100644 (file)
--- a/MolLib/Dcd.h
+++ b/MolLib/Dcd.h
@@ -32,6 +32,7 @@ extern "C" {
                
 typedef int Int32;
 typedef float SFloat32;
+typedef double SFloat64;
 
 typedef struct DcdRecord {
     int fd;                 /*  File descripter  */    
@@ -41,10 +42,14 @@ typedef struct DcdRecord {
        Int32 nstart;           /*  Start timestep  */
        Int32 ninterval;        /*  Number of timesteps between frames  */
        Int32 nend;             /*  Last timestep  */
-       Int32 with_unitcell;    /*  Has a unit cell information?  */
+       Int32 nfix;             /*  Number of fixed atoms  */
+       Int32 n4dim;            /*  Charmm 4th dimension  */
+       Int32 ncharmver;        /*  Charmm version  */
+       Int32 nextra;           /*  Has Charmm extra block  */
     SFloat32 delta;         /*  Step time  */
     SFloat32 globalcell[6]; /*  cell size and origin; used when with_unitcell == 0 */
     off_t header_size;      /*  Header size  */
+       off_t block_size;       /*  Block size (set after reading/writing the first frame)  */
 } DcdRecord;
 
 int DcdOpen(const char *name, DcdRecord *dr);

diff --git a/MolLib/Molecule.c b/MolLib/Molecule.c
index 016c4a2..8e12129 100755 (executable)
--- a/MolLib/Molecule.c
+++ b/MolLib/Molecule.c
@@ -3082,7 +3082,7 @@ MoleculeReadCoordinatesFromDcdFile(Molecule *mp, const char *fname, char *errbuf
 {
        DcdRecord dcd;
        SFloat32 *xp, *yp, *zp;
-       Vector *vp;
+       Vector *vp, *cp;
        IntGroup *ig;
        int n;
        errbuf[0] = 0;
@@ -3115,6 +3115,9 @@ MoleculeReadCoordinatesFromDcdFile(Molecule *mp, const char *fname, char *errbuf
        }
 
        vp = (Vector *)calloc(sizeof(Vector), mp->natoms * dcd.nframes);
+       if (dcd.nextra)
+               cp = (Vector *)calloc(sizeof(Vector), dcd.nframes * 4);
+       else cp = NULL;
        xp = (SFloat32 *)malloc(sizeof(SFloat32) * dcd.natoms);
        yp = (SFloat32 *)malloc(sizeof(SFloat32) * dcd.natoms);
        zp = (SFloat32 *)malloc(sizeof(SFloat32) * dcd.natoms);
@@ -3122,6 +3125,7 @@ MoleculeReadCoordinatesFromDcdFile(Molecule *mp, const char *fname, char *errbuf
        if (vp == NULL || xp == NULL || yp == NULL || zp == NULL || ig == NULL) {
                snprintf(errbuf, errbufsize, "Cannot allocate memory");
                if (vp) free(vp);
+               if (cp) free(cp);
                if (xp) free(xp);
                if (yp) free(yp);
                if (zp) free(zp);
@@ -3131,7 +3135,8 @@ MoleculeReadCoordinatesFromDcdFile(Molecule *mp, const char *fname, char *errbuf
        for (n = 0; n < dcd.nframes; n++) {
                int i;
                Vector *vpp;
-               if (DcdReadFrame(&dcd, n, xp, yp, zp, dcd.globalcell)) {
+               SFloat32 dcdcell[6];
+               if (DcdReadFrame(&dcd, n, xp, yp, zp, dcdcell)) {
                        snprintf(errbuf, errbufsize, "Read error in dcd file");
                        goto exit;
                }
@@ -3140,29 +3145,44 @@ MoleculeReadCoordinatesFromDcdFile(Molecule *mp, const char *fname, char *errbuf
                        vpp->y = yp[i];
                        vpp->z = zp[i];
                }
-       }
-       /*  TODO: implement frame-specific cells  */
-       if (MoleculeInsertFrames(mp, ig, vp, NULL) < 0)
+               if (cp != NULL) {
+                       Double sing;
+                       vpp = &cp[n * 4];
+                       /*  dcdcell = {a, gamma, b, beta, alpha, c} */
+                       /*  angles are described either in cosines (Charmm and NAMD > 2.5) or degrees (NAMD 2.5)  */
+                       if (dcdcell[1] < -1.0 || dcdcell[1] > 1.0 || dcdcell[3] < -1.0 || dcdcell[3] > 1.0 || dcdcell[4] < -1.0 || dcdcell[4] > 1.0) {
+                               dcdcell[4] = cos(dcdcell[4] * kDeg2Rad);  /*  cos(alpha)  */
+                               dcdcell[3] = cos(dcdcell[3] * kDeg2Rad);  /*  cos(beta)  */
+                               dcdcell[1] = cos(dcdcell[1] * kDeg2Rad);  /*  cos(gamma)  */
+                       }
+                       /*  a axis lies along the cartesian x axis  */
+                       sing = sqrt(1 - dcdcell[1] * dcdcell[1]);
+                       vpp[0].x = dcdcell[0];
+                       vpp[0].y = 0;
+                       vpp[0].z = 0;
+                       vpp[1].x = dcdcell[2] * dcdcell[1];
+                       vpp[1].y = dcdcell[2] * sing;
+                       vpp[1].z = 0;
+                       vpp[2].x = dcdcell[5] * dcdcell[3];
+                       vpp[2].y = dcdcell[5] * (dcdcell[4] - dcdcell[3] * dcdcell[1]) / sing;
+                       vpp[2].z = sqrt(dcdcell[5] * dcdcell[5] - vpp[2].x * vpp[2].x - vpp[2].y * vpp[2].y);
+                       vpp[3].x = vpp[3].y = vpp[3].z = 0.0;
+                       if (mp->cell == NULL) {
+                               /*  Create periodicity if not present  */
+                               static const char pflags[] = {1, 1, 1};
+                               MoleculeSetPeriodicBox(mp, &vpp[0], &vpp[1], &vpp[2], &vpp[3], pflags);
+                       }
+               }
+       }
+       if (MoleculeInsertFrames(mp, ig, vp, cp) < 0)
                snprintf(errbuf, errbufsize, "Cannot insert frames");
-       if (dcd.with_unitcell) {
-               Vector ax, ay, az, orig;
-               char flags[3] = {1, 1, 1};
-               ax.x = dcd.globalcell[0]; ax.y = ax.z = 0;
-               ay.y = dcd.globalcell[1]; ay.x = ay.z = 0;
-               az.z = dcd.globalcell[2]; az.x = az.y = 0;
-               orig.x = dcd.globalcell[3];
-               orig.y = dcd.globalcell[4];
-               orig.z = dcd.globalcell[5];
-               if (MoleculeSetPeriodicBox(mp, &ax, &ay, &az, &orig, flags) != 0) {
-                       snprintf(errbuf, errbufsize, "Cannot set unit cell");
-                       goto exit;
-               }
-       }
        mp->startStep = dcd.nstart;
        mp->stepsPerFrame = dcd.ninterval;
        mp->psPerStep = dcd.delta;
 exit:
        DcdClose(&dcd);
+       if (cp != NULL)
+               free(cp);
        free(vp);
        free(xp);
        free(yp);
@@ -3753,18 +3773,12 @@ MoleculeWriteToDcdFile(Molecule *mp, const char *fname, char *errbuf, int errbuf
        if (dcd.ninterval == 0)
                dcd.ninterval = 1;
        dcd.nend = dcd.nstart + (dcd.nframes - 1) * dcd.ninterval;
-       if (mp->cell != NULL) {
-               dcd.with_unitcell = 1;
-               dcd.globalcell[0] = VecLength(mp->cell->axes[0]);
-               dcd.globalcell[1] = VecLength(mp->cell->axes[1]);
-               dcd.globalcell[2] = VecLength(mp->cell->axes[2]);
-               dcd.globalcell[3] = mp->cell->origin.x;
-               dcd.globalcell[4] = mp->cell->origin.y;
-               dcd.globalcell[5] = mp->cell->origin.z;
-       }
+       if (mp->cell != NULL)
+               dcd.nextra = 1;
        dcd.delta = mp->psPerStep;
        if (dcd.delta == 0.0)
                dcd.delta = 1.0;
+       dcd.ncharmver = 24;
        n = DcdCreate(fname, &dcd);
        if (n != 0) {
                if (n < 0)
@@ -3805,6 +3819,15 @@ MoleculeWriteToDcdFile(Molecule *mp, const char *fname, char *errbuf, int errbuf
                        memset(yp + i, 0, sz);
                        memset(zp + i, 0, sz);
                }
+               if (n < mp->nframes && mp->frame_cells != NULL) {
+                       Vector *cp = &(mp->frame_cells[n * 4]);
+                       dcd.globalcell[0] = VecLength(cp[0]);
+                       dcd.globalcell[2] = VecLength(cp[1]);
+                       dcd.globalcell[5] = VecLength(cp[2]);
+                       dcd.globalcell[1] = VecDot(cp[0], cp[1]) / (dcd.globalcell[0] * dcd.globalcell[2]);
+                       dcd.globalcell[3] = VecDot(cp[0], cp[2]) / (dcd.globalcell[0] * dcd.globalcell[5]);
+                       dcd.globalcell[4] = VecDot(cp[1], cp[2]) / (dcd.globalcell[2] * dcd.globalcell[5]);                     
+               }                       
                if (DcdWriteFrame(&dcd, n, xp, yp, zp, dcd.globalcell)) {
                        snprintf(errbuf, errbufsize, "Write error in dcd file");
                        goto exit;

diff --git a/wxSources/MyDocument.cpp b/wxSources/MyDocument.cpp
index 7ebb2e1..6e5b53a 100755 (executable)
--- a/wxSources/MyDocument.cpp
+++ b/wxSources/MyDocument.cpp
@@ -273,6 +273,7 @@ MyDocument::OnImport(wxCommandEvent& event)
                }
                /*  Insert Import-only file types before "All files"  */
                wildcard += _T("|AMBER mdcrd file (*.crd;*.mdcrd)|*.crd;*.mdcrd");
+               wildcard += _T("|DCD file (*.dcd)|*.dcd");
                if (i == -1)
                        wildcard += (_T("|") + desc + _T(" (") + filter + _T(")|") + filter);
        }
@@ -294,22 +295,35 @@ void
 MyDocument::OnExport(wxCommandEvent& event)
 {
        wxString wildcard;
+       wxFileName fname(GetFilename());
+       wxString fnstr;
+       GetPrintableName(fnstr);
        {
                /*  File filter is built from MyDocManager information  */
                wxString desc, filter, ext;
                int i;
                MyDocManager *docm = wxGetApp().DocManager();
+               if ((i = fnstr.Find('.', true)) != wxNOT_FOUND) {
+                       fnstr = fnstr.Mid(0, i);
+               }
                for (i = 0; docm->GetDocumentDescriptionAtIndex(i, &desc, &filter, &ext); i++) {
-                       if (ext == _T("out") || ext == _T("log") || ext == _T("fchk"))
+                       if (ext == _T("mbsf") || ext == _T("out") || ext == _T("log") || ext == _T("fchk"))
                                continue;
+                       if (filter.Contains(_T("*.*"))) {
+                               i = -1;
+                               break;
+                       }
                        if (wildcard != _T("")) {
                                wildcard += (_T("|"));
                        }
                        wildcard += (desc + _T(" (") + filter + _T(")|") + filter);
                }
+               wildcard += _T("|AMBER mdcrd file (*.crd;*.mdcrd)|*.crd;*.mdcrd");
+               wildcard += _T("|DCD file (*.dcd)|*.dcd");
+               if (i == -1)
+                       wildcard += (_T("|") + desc + _T(" (") + filter + _T(")|") + filter);
        }
-       
-       wxFileDialog *dialog = new wxFileDialog(NULL, _T(""), _T(""), _T(""), wildcard, wxFD_SAVE | wxFD_OVERWRITE_PROMPT |   wxFD_CHANGE_DIR);
+       wxFileDialog *dialog = new wxFileDialog(NULL, _T("Export coordinates"), fname.GetPath(), fnstr + _T(".psf"), wildcard, wxFD_SAVE | wxFD_OVERWRITE_PROMPT | wxFD_CHANGE_DIR);
        if (dialog->ShowModal() == wxID_OK) {
                char *p = strdup((const char *)(dialog->GetPath().mb_str(wxConvFile)));
                MoleculeLock(mol);