Formatting the HISTORY files of DL_MESO_DPD

Purpose of Module

This module format_history.f90 is a post-processing utility for DL_MESO_DPD, the Dissipative Particle Dynamics (DPD) code from the DL_MESO package.

It converts the trajectory (HISTORY) file from unformatted to a human readable form, (optionally) including explanatory comments about all the quantities. This module is mainly for learning/checking purposes. The first aim is to help the user to check that the system was prepared as intended (e.g., showing all the bead properties and initial positions, all the bonds etc.). The idea is to use it on small systems when familiarizing with the structure of input files needed for the simulation. Secondly, it can be used as a starting point for user-defined analyses of trajectories.

Background Information

The base code for this module is DL_MESO_DPD, the Dissipative Particle Dynamics code from the mesoscopic simulation package DL_MESO, developed by M. Seaton at Daresbury Laboratory. This open source code is available from STFC under both academic (free) and commercial (paid) licenses. The module is to be used with DL_MESO in its currently released version, version 2.7 (dating December 2018).

A version of the module that works with HISTORY files generated by the previous version of DL_MESO, version 2.6 (dating November 2015), can be found in the old-v2.6 directory.

Testing

The present module is compiled with the available Fortran 2003 compiler, e.g.:

gfortran -o format.exe format_history.f90

and the executable must be in the same directory of the HISTORY file to be analyzed. To test the module, run the simulation with the toy input files given in the following. (Note that these files contain commented lines as suggestions for further tests.) For the CONTROL file

Simple test

volume 3.0  3.0  3.0
temperature 1.0
cutoff 1.0

timestep 0.01
steps 6
equilibration steps 2
traj  2  2  0
stats every 2
stack size 2
print every 2
job time 100.0
close time 10.0

#surface shear y
#surface frozen  x
#surface hard  x

ensemble nvt mdvv

finish

and for the FIELD file

Simple test

SPECIES 3
A    1.0  0.0   1   0
B    1.0  0.0   0   0
C    1.0  0.0   0   0 

MOLECULES 2
AB
nummols 1
beads 2
A    0.0 0.0 0.0
B    0.1 0.0 0.0
bonds 1
harm 1 2 5.0 0.0
finish
AC
nummols 1
beads 2
A    0.0 0.0 0.0
C    0.1 0.0 0.0
bonds 1
harm 1 2 3.0 0.0
finish

INTERACTIONS 3
A    A    dpd 25.0 1.0 4.5
B    B    dpd 25.0 1.0 4.5
C    C    dpd 25.0 1.0 4.5

#EXTERNAL
#shear  3.0  0.0  0.0

CLOSE

After analyzing the trajectories, for a serial run (i.e., on a single processor core) and for lcomm , lmcheck and sorted all set to .TRUE., this output should be printed on the screen

 # Check of beads: i, ltp(i), ltm(i), mole(i)
           1           1           0           0
           2           1           1           1
           3           2           1           1
           4           1           2           2
           5           3           2           2
 # Check of molecules: nammol(i), nbdmol(i), nbomol(i), nmol(i)
 AB                 2           1           1
 AC                 2           1           1
 # Total number of molecules =            2
 # Check of bonds: bndbtl(i,1), bndbtl(i,2)
           2           3
           4           5

and the HISTORY-F file should be

 # Simulation name:
 Simple test                                                                     
 # nspe, nmoldef, nusyst, nsyst, numbond
           3           2           1           5           2
 # keytrj, srfx, srfy, srfz
           0           0           0           0
 # SPECIES:
 # namspe, amass, rcii, chge, lfrzn
 A               1.000        1.000        0.000    0
 B               1.000        1.000        0.000    0
 C               1.000        1.000        0.000    0
 # MOLECULES:
 # nammol
 AB      
 AC      
 # BEADS:
 # global, species, molecule, chain
           1           1           0           0
           2           1           1           1
           3           2           1           1
           4           1           2           2
           5           3           2           2
 # BONDS:
 # extremes of the bond
           2           3
           4           5
 # --- TRAJECTORIES --- (key =           0 )
 # mglobal, x, y, z
 # time, nbeads, dimx, dimy, dimz, shrdx, shrdy, shrdz
     0.000             5     3.000        3.000        3.000        0.000        0.000        0.000
 # snapshot number:           1
         1  -5.594619E-03   -1.752832E-03    6.889737E-03
         2  -1.691781E-01    6.786290E-01    1.166801E-02
         3  -2.345695E-01    7.399172E-01   -1.447126E-01
         4   5.198776E-01   -8.222446E-01    1.038467E-02
         5   4.631501E-01   -8.058678E-01    1.464801E-01
 # time, nbeads, dimx, dimy, dimz, shrdx, shrdy, shrdz
     0.020             5     3.000        3.000        3.000        0.000        0.000        0.000
 # snapshot number:           2
         1  -1.061544E-02   -6.284880E-03    1.390368E-02
         2  -1.395560E-01    6.796778E-01    4.361174E-02
         3  -2.443224E-01    7.728033E-01   -1.800995E-01
         4   5.280036E-01   -8.191712E-01   -3.253643E-02
         5   4.401757E-01   -8.383441E-01    1.858304E-01
 # time, nbeads, dimx, dimy, dimz, shrdx, shrdy, shrdz
     0.040             5     3.000        3.000        3.000        0.000        0.000        0.000
 # snapshot number:           3
         1  -1.483900E-02   -1.635085E-02    2.109335E-02
         2  -1.083779E-01    6.830112E-01    8.086194E-02
         3  -2.571191E-01    8.101380E-01   -2.211317E-01
         4   5.390852E-01   -8.157553E-01   -8.230434E-02
         5   4.149363E-01   -8.723621E-01    2.321906E-01

If the simulation is run in parallel, the particles may not necessarily be written to the HISTORY file in order of particle index, but the module can sort the particles in each trajectory snapshot before printing to the HISTORY-F file.

Source Code

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
PROGRAM format_history
!***********************************************************************************
!
! module to format dl_meso HISTORY files
!
! authors - m. a. seaton & s. chiacchiera, february 2017 (amended january 2021)
!
!**********************************************************************************
      IMPLICIT none
      INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND (15, 307)
      INTEGER, PARAMETER :: li = SELECTED_INT_KIND (12)
      INTEGER, PARAMETER :: ntraj=10,nuser=5
      INTEGER, PARAMETER :: endversion = 1

      CHARACTER(80) :: text
      CHARACTER(8), ALLOCATABLE :: namspe (:), nammol (:)

      INTEGER, ALLOCATABLE :: ltp (:), ltm (:), mole (:), bndtbl (:,:)
      INTEGER, ALLOCATABLE :: nbdmol (:), nbomol (:), readint (:), globindex (:)
      INTEGER :: chain, imol, ioerror, i, k, nmoldef, numframe
      INTEGER :: nspe, nbeads, nusyst, nsyst, global, species, molecule, numbond
      INTEGER :: nummol, lfrzn, rnmol, keytrj, srfx, srfy, srfz
      INTEGER :: bead1, bead2
      INTEGER :: nform
      INTEGER :: endver, Dlen, nstep, framesize, lend, leni
      INTEGER(KIND=li) :: filesize, mypos, headerpos, currentpos, lend_li, leni_li, framesizeli, numbeadsli

      REAL(KIND=dp), ALLOCATABLE :: nmol (:), readdata (:)
      REAL(KIND=dp) :: dimx, dimy, dimz, shrdx, shrdy, shrdz
      REAL(KIND=dp) :: amass, rcii, chge
      REAL(KIND=dp) :: time

      LOGICAL :: eof, lcomm, lmcheck, swapend, bigend, sorted

      ! Switches for commenting, checking molecules and sorting particles in output
      
      lcomm =   .TRUE.
      lmcheck = .TRUE.
      sorted  = .TRUE.

!     determine number of bytes for selected double precision kind
!     (the default SELECTED_REAL_KIND (15, 307) should return 8 bytes)

      lend = STORAGE_SIZE (1.0_dp) / 8
      leni = BIT_SIZE (1) / 8
      lend_li = INT (lend, KIND=li)
      leni_li = INT (leni, KIND=li)

!     check endianness of machine

      bigend = (IACHAR(TRANSFER(1,"a"))==0)

      ! Determine if HISTORY file exists, which endianness to use,
      ! if type of real is correct

      INQUIRE (file = 'HISTORY', EXIST = eof)
      IF (.NOT. eof) THEN
        PRINT *, "ERROR: cannot find HISTORY file"
        STOP
      END IF

      OPEN (ntraj, file = 'HISTORY', access = 'stream', form = 'unformatted', status = 'unknown')

      swapend = .false.
      READ (ntraj) endver, Dlen

      IF (endver/=endversion) THEN
        swapend = .true.
        CLOSE (ntraj)
        IF (bigend) THEN
          OPEN (ntraj, file = 'HISTORY', access = 'stream', form = 'unformatted', status = 'unknown', convert = 'little_endian')
        ELSE
          OPEN (ntraj, file = 'HISTORY', access = 'stream', form = 'unformatted', status = 'unknown', convert = 'big_endian')
        END IF
        READ (ntraj) endver, Dlen
        IF (endver/=endversion) THEN
          PRINT *, "ERROR: corrupted HISTORY file or created with incorrect version of DL_MESO"
          STOP
        END IF
      END IF

      IF (Dlen/=lend) THEN
        PRINT *, "ERROR: incorrect type of real number used in HISTORY file"
        PRINT *, "       recompile format_history.f90 with reals of ", Dlen, " bytes"
        STOP
      END IF

      ! Open the output file
      nform = ntraj + 1
      OPEN (nform, file = 'HISTORY'//"-F", status = 'replace')

!     read file size, number of frames and timestep numbers

      READ (ntraj) filesize, numframe, nstep

      ! read the number of beads, molecules and bonds
      ! Arrays are filled with names of particles and molecules: if checking molecules,
      ! arrays for species, molecule types etc. also filled

      READ (ntraj) text
      READ (ntraj) nspe, nmoldef, nusyst, nsyst, numbond, keytrj, srfx, srfy, srfz

      IF (lcomm) WRITE (nform,*) "# Simulation name:"
      WRITE (nform,*) text

      IF (lcomm) WRITE (nform,*) "# nspe, nmoldef, nusyst, nsyst, numbond"
      WRITE (nform,*) nspe, nmoldef, nusyst, nsyst, numbond
      IF (lcomm) WRITE (nform,*) "# keytrj, srfx, srfy, srfz"
      WRITE (nform,*) keytrj, srfx, srfy, srfz

      framesize = (keytrj+1) * 3
      ALLOCATE (namspe (nspe), nammol (nmoldef), globindex (nsyst), readint (nsyst), readdata (framesize))
      IF (lmcheck) THEN 
         ALLOCATE (ltp (1:nsyst), ltm (1:nsyst), mole (1:nsyst))
         ALLOCATE (nmol (1:nmoldef), nbdmol (1:nmoldef), nbomol (1:nmoldef))
         ALLOCATE (bndtbl (numbond, 2))
      END IF
      
      IF (lcomm) WRITE (nform,*) "# SPECIES:"
      IF (lcomm) WRITE (nform,*) "# namspe, amass, rcii, chge, lfrzn"
      DO i = 1, nspe
        READ (ntraj) namspe (i), amass, rcii, chge, lfrzn
        WRITE (nform,96) namspe (i), amass, rcii, chge, lfrzn
      END DO

      IF (nmoldef>0) THEN
        IF (lcomm) WRITE (nform,*) "# MOLECULES:"
        IF (lcomm) WRITE (nform,*) "# nammol"
        DO i = 1, nmoldef
          READ (ntraj) nammol (i)
          WRITE (nform,*) nammol (i)
        END DO
      END IF

      ! (if required) read and fill arrays with properties of beads and molecules

      nummol = 0 ! counter for number of molecules

      IF (lcomm) WRITE (nform,*) "# BEADS:"
      IF (lcomm) WRITE (nform,*) "# global, species, molecule, chain"
      IF (lmcheck) THEN
        ! Build ltp, ltm, mole
        DO i = 1, nsyst
          READ (ntraj) global, species, molecule, chain
          ltp (global) = species
          ltm (global) = molecule
          mole (global) = chain
          nummol = MAX (nummol, chain)
          WRITE (nform,*) global, species, molecule, chain
        END DO
      ELSE
        DO i = 1, nsyst
          READ (ntraj) global, species, molecule, chain
          WRITE (nform,*) global, species, molecule, chain
        END DO
      END IF

      IF (numbond>0) THEN
        IF (lcomm) WRITE (nform,*) "# BONDS:"
        IF (lcomm) WRITE (nform,*) "# extremes of the bond"
        IF (lmcheck) THEN
         ! Build bndtbl
          DO i = 1, numbond
            READ (ntraj) bead1, bead2
            bndtbl (i, 1) = bead1
            bndtbl (i, 2) = bead2
            WRITE (nform,*) bead1, bead2
          END DO
        ELSE
          DO i = 1, numbond
            READ (ntraj) bead1, bead2
            WRITE (nform,*) bead1, bead2
          END DO
        END IF
      END IF

!     reached end of header: find current position in file

      INQUIRE (unit=ntraj, POS=headerpos)
      framesizeli = INT (framesize, KIND=li)
      numbeadsli = INT (nsyst, KIND=li)

      IF (lmcheck) THEN
      ! determine numbers of molecules, beads and bonds per molecule type
         nmol = 0.0_dp
         nbdmol = 0
         nbomol = 0 
         chain = 0
         imol = 0 ! necessary to avoid out of bounds
         
         DO i = 1, nsyst
            IF (mole (i) /= chain) THEN
               chain = mole (i)
               imol = ltm (i)
               nmol (imol) = nmol (imol) + 1.0_dp
            END IF
            IF (imol > 0) nbdmol (imol) = nbdmol (imol) + 1
         END DO

         DO i = 1, numbond
            imol = ltm (bndtbl (i,1))
            nbomol (imol) = nbomol (imol) + 1
         END DO
            
         DO i = 1, nmoldef
            rnmol = NINT (nmol (i))
            IF (rnmol>0) THEN
               nbdmol (i) = nbdmol (i) / rnmol
               nbomol (i) = nbomol (i) / rnmol
            END IF
         END DO

         ! Write to std output the arrays built
         WRITE (*,*) "# Check of beads: i, ltp(i), ltm(i), mole(i)"
         DO i = 1, nsyst
            WRITE(*,*) i, ltp (i), ltm (i), mole (i)
         END DO

         !Check of molecule beads and numbers
         IF (nmoldef>0) THEN
            WRITE (*,*) "# Check of molecules: nammol(i), nbdmol(i), nbomol(i), nmol(i)"
            DO i = 1, nmoldef
               WRITE (*,*) nammol (i), nbdmol (i), nbomol (i), NINT(nmol(i))
            END DO
            WRITE (*,*) "# Total number of molecules = ",nummol
         END IF

         ! Write to std output bndtbl
         IF (numbond > 0) THEN
            WRITE (*,*) "# Check of bonds: bndbtl(i,1), bndbtl(i,2)"
            DO i = 1, numbond
               WRITE (*,*) bndtbl (i,1), bndtbl (i,2)
            END DO
         END IF
      END IF

      ! Now read in trajectories

      eof = .false.

      IF (lcomm) WRITE (nform,*) "# --- TRAJECTORIES --- (key =", keytrj,")"
      SELECT CASE (keytrj)
      CASE (0)
        IF (lcomm) WRITE (nform,*) "# mglobal, x, y, z"
      CASE(1)
        IF (lcomm) WRITE (nform,*) "# mglobal, x, y, z, vx, vy, vz"
      CASE(2)
        IF (lcomm) WRITE (nform,*) "# mglobal, x, y, z, vx, vy, vz, fx, fy, fz"
      END SELECT

      DO k = 1, numframe

        currentpos = headerpos + INT (k-1, KIND=li) * ((7_li + numbeadsli * framesizeli) * lend_li + (1_li + numbeadsli) * leni_li)
        READ (ntraj, POS=currentpos, IOSTAT=ioerror) time, nbeads, dimx, dimy, dimz, shrdx, shrdy, shrdz

        IF (ioerror/=0) THEN
          eof = .true.
          IF (k==1) THEN
            PRINT *, 'ERROR: cannot find trajectory data in HISTORY file'
            STOP
          END IF
          EXIT
        END IF

        IF (lcomm) WRITE (nform,*) "# time, nbeads, dimx, dimy, dimz, shrdx, shrdy, shrdz"
        WRITE (nform,98) time, nbeads, dimx, dimy, dimz, shrdx, shrdy, shrdz
            
        IF (lcomm) WRITE (nform,*) "# snapshot number:", k

        IF (sorted) THEN

          READ (ntraj) readint (1:nbeads)
          CALL quicksort_integer_indexed (readint, 1, nbeads, globindex)
          DO i = 1, nbeads
            global = globindex (i)
            mypos = currentpos + leni_li * (1_li + numbeadsli) + (7_li + INT (global-1, KIND=li) * framesizeli) * lend_li
            READ (ntraj, POS=mypos) readdata (1:framesize)
            WRITE (nform,99) global, readdata (1:framesize)
          END DO

        ELSE

          READ (ntraj) globindex (1:nbeads)
          DO i = 1, nbeads
            READ (ntraj, POS=mypos) readdata (1:framesize)
            WRITE (nform,99) global, readdata (1:framesize)
          END DO

        END IF

      END DO
         
      ! Close the trajectory file
      CLOSE (ntraj)

      ! close the output file
      CLOSE (nform)

      DEALLOCATE (readint, readdata, globindex)
      DEALLOCATE (namspe, nammol)
      IF (lmcheck) DEALLOCATE (ltp, ltm, mole, nmol, nbdmol, bndtbl, nbomol)

99    FORMAT(I10,2x,1p,9(e13.6,3x))
98    FORMAT(f10.3,3x,1x,I10,6(f10.3,3x))
96    FORMAT(A9,3x,3(f10.3,3x),I2)

CONTAINS

      SUBROUTINE quicksort_integer_indexed (list, stride, n, indices)

!**********************************************************************
!
!     sort integers in array into numerical order, recording original
!     positions of values (routine to prepare indices array)
!
!     copyright ukri stfc daresbury laboratory
!     authors - m. a. seaton august 2013
!
!**********************************************************************

      INTEGER, INTENT (INOUT) :: list (:)
      INTEGER, INTENT (IN) :: stride, n
      INTEGER, INTENT (OUT) :: indices (:)
      INTEGER :: i

      DO i = 1, n
        indices (i) = i
      END DO

      CALL qsort_integer (list, indices, stride, 1, n)

      END SUBROUTINE quicksort_integer_indexed

      RECURSIVE SUBROUTINE qsort_integer (list, index, stride, low, high)

!**********************************************************************
!
!     sort integers in array into numerical order, recording original
!     positions of values
!
!     copyright ukri stfc daresbury laboratory
!     authors - m. a. seaton august 2013
!
!**********************************************************************

      INTEGER, INTENT (INOUT) :: list (:), index (:)
      INTEGER, INTENT (IN) :: low, high
      INTEGER, INTENT (IN) :: stride
      INTEGER :: i, j, k, reference, temp

      IF (high < low + 6) THEN

!     resort to bubble sort for very small lists (5 items or fewer)

        DO i = low, high - 1
          DO j = i+1, high
            IF (list (stride * (i - 1) + 1) > list (stride * (j - 1) + 1)) THEN
              DO k = 1, stride
                temp = list (stride * (i - 1) + k)
                list (stride * (i - 1) + k) = list (stride * (j - 1) + k)
                list (stride * (j - 1) + k) = temp
              END DO
              temp = index (i)
              index (i) = index (j)
              index (j) = temp
            END IF
          END DO
        END DO

      ELSE

!     apply partition-based sort

        reference = list (stride * ((low+high)/2 - 1) + 1)
        i = low - 1
        j = high + 1
        DO
          DO
            i = i + 1
            IF (list (stride * (i-1) + 1) >= reference) EXIT
          END DO
          DO
            j = j - 1
            IF (list (stride * (j-1) + 1) <= reference) EXIT
          END DO
          IF (i < j) THEN
            DO k = 1, stride
              temp = list (stride * (i-1) + k)
              list (stride * (i-1) + k) = list (stride * (j-1) + k)
              list (stride * (j-1) + k) = temp
            END DO
            temp = index (i)
            index (i) = index (j)
            index (j) = temp
          ELSE IF (i == j) THEN
            i = i + 1
            EXIT
          ELSE
            EXIT
          END IF
        END DO

        IF (low<j) CALL qsort_integer (list, index, stride, low, j)
        IF (i<high) CALL qsort_integer (list, index, stride, i, high)

      END IF

      END SUBROUTINE qsort_integer

END PROGRAM format_history