Dear Colleague,
Five bugs have been spotted in the Dissipative Particle Dynamics (DPD) code of
DL_MESO version 2.6.
The application of Ewald sums for Slater-based charge smearing is currently
incorrect: as the code currently stands, the real space contributions to
potentials and forces are higher than they should be for shorter distances
between ion pairs. The original (incorrect) formulation of the real-space terms
has also led to incorrect terms for potential and force corrections on pairs of
charged frozen particles. (Many thanks are due to Patrick Warren at Unilever for
bringing this bug to our attention.) Long-range effects due to the reciprocal
space part of the Ewald sum are not affected by this bug and the real space
potentials and forces for Gaussian charge smearing have been applied correctly.
To fix this bug, changes need to be made to ewald_module.f90 and
ewald_module_omp.f90: line numbers for the latter are shown in brackets. In the
subroutine ewald_real_slater in ewald_module.f90 (ewald_module_omp.f90), line
217 (242) needs to be changed to:
scrn = EXP (-2.0_dp * brr) / rrr
line 219 (244) to:
epot = bjerelec * qiqj * (erfr - (1.0_dp + brr) * scrn)
and lines 223 and 224 (247 and 248) to:
eforce = bjerelec * qiqj * (2.0_dp * alphaew / rtpi * EXP (-arr * arr) + erfr - &
&scrn * (1.0_dp + 2.0_dp * brr * (1.0_dp + brr))) / rsq
The subroutine ewald_frozen_slater includes three sets of calculations for the
correction terms due to frozen particle pairs, and all three need to be
corrected. Lines 931, 932, 963, 964, 1088 and 1089 in ewald_module.f90 (lines
1063, 1064, 1095, 1096, 1140 and 1141 in ewald_module_omp.f90 can all be
commented out as the variables brr and scrn are no longer needed, lines 934, 966
and 1011 (1066, 1098 and 1143) each need to be changed to:
epot = bjerelec * qiqj * erfr
and lines 936-937, 968-969 and 1013-1014 (1068-1069, 1100-1101 and 1145-1146)
should be replaced with:
eforce = bjerelec * qiqj * (erfr - 2.0_dp * alphaew / rtpi * EXP (-arr * arr)) / rsq
Additionally, in ewald_module_omp.f90, lines 170 and 400 can be changed to the
following to improve efficiency:
!$OMP& bjerelec,idnode,relec,rel2,lfrzn,chge,tnum,xxyyzz,mxpcellew) &
The user manual has also been modified to reflect the change in application of
the Ewald sum for Slater-based smearing, and the demonstration case
Polyelectrolyte has been updated to use the correct potentials and forces.
When calculations in serial are being carried out without many-body DPD
interactions, the code may encounter a segmentation fault due to erroneous
determination of the numbers of many-body DPD link cells: this error does not
occur when running the parallel version of the code. (Many thanks are due to
Neil Morgan at STFC Hartree Centre for drawing our attention to this bug.) To
fix this bug, change line 201 in domain_module_ser.f90 to:
IF (lmb .AND. nlmbx2*nlmby2*nlmbz2/=nlold) THEN
When carrying out calculations in parallel involving a barostat (e.g. using
constant pressure ensembles), if the initial system pressure is considerably
different to the required value, the simulation may halt with an error message
indicating the communication buffer arrays are too small. This bug can be fixed
by changing line 4099 in domain_module.f90 to:
IF (lmb .OR. btype>0 .OR. step