Dear Colleague,

A number of bugs in both the Lattice Boltzmann and Dissipative Particle Dynamics 
codes have been discovered.

One bug was discovered in the Lattice Boltzmann code: multiple phase systems 
with zero densities for particular species at certain lattice points may cause 
division-by-zero issues when applying velocity corrections for external 
(interfacial) forces for BGK collisions. (Many thanks are due to Ernesto Monaco 
for bringing this bug to our attention.)

Changing line 1125 in lbpSUB.cpp (in the routine fSiteCollisionBGK) to:

      speed[i] = sitespeed[i];
      if(onemass>0.0) speed[i] += (*pt3) / (onemass * lbtf[j]);

and line 1224 in the same file (fSiteCollisionBGKGuo) to:

      speed[i] = sitespeed[i];
      if(onemass>0.0) speed[i] += 0.5 * force[i] / onemass;

will fix this bug. Note that this bug should never occur for incompressible BGK 
or multiple-relaxation-time collisions.

Four bugs were discovered in the Dissipative Particle Dynamics code:

1. The routines contracting bond, angle and dihedral tables are unnecessarily 
   zeroing the tables after the remaining bonds (angles, dihedrals) are shifted 
   to the beginning. This bug has no effect on the accuracy of calculations but 
   has been observed to cause computational inefficiencies, particularly when 
   running with large numbers of processors and/or polymers in the system. 
   Commenting out or removing lines 274-282, 317-326 and 364-372 in 
   bond_module.f90 will fix this bug. (Note that the lines reassigning numbers 
   of bonds, angles and dihedrals MUST be retained.)

2. The readconfig routine in start_module.f90 currently assumes the molecules 
   are to be read in sequential order (i.e. it cannot assign molecule numbers 
   correctly if the global bead numbers are not in sequential order) and will 
   not assign dihedral data correctly. The first of these problems can be fixed 
   by changing line 996 in start_module.f90 to:

        imol = 0
        ntop = 0
        j = 0

        DO WHILE (global>molstart(j+1))
          j = j + 1
        END DO

   while the second can be fixed by changing line 1141 to:

        global = molstart (j) + dhdinp2 (imol, k)

   Lines 961 to 963 can also be commented out or deleted.

3. When using global bond interactions in parallel, only the HISTORY* file for 
   processor 0 is assigned bond data, but the total number of bonds is given on 
   all files. This causes the utilities that read HISTORY* files to attempt to 
   read bond data that does not exist. To fix this problem, line 659 in 
   start_module.f90 should be changed to:

        IF (lgbnd .AND. idnode>0) THEN
          WRITE (ntraj) isbsv, nspe, nmoldef, nusyst, nsyst, nbeads, 0
        ELSE
          WRITE (ntraj) isbsv, nspe, nmoldef, nusyst, nsyst, nbeads, nbonds
        END IF

4. Both the traject.f90 and trajectselected.f90 utilities may produce unreadable 
   traject.vtf files for systems without bonds if the final numbered bead is of 
   the most common type. To fix this bug, line 142 in traject.f90 and line 153 
   in trajectselected.f90 should be changed to:

        IF ((ltp (i)/=imxspe .OR. i==last) .AND. ltm(i)==0) THEN

   which will force the utilities to write data for the last numbered bead.

We strongly recommend that all registered users should either make the code 
changes listed above or download DL_MESO again using the instructions they 
received on registration to obtain the corrected version of the software. Re-
registration is not required and the same password is used for decryption.

Michael Seaton
2011-07-28