Bond/angle style Gaussian

tcnicholas · November 18, 2022, 10:14am

Hi. I am using the Gaussian bond and angle styles in a simulation. Comparing the expression in the LAMMPS documentation to the original paper (J. Polym. Sci. B Polym. Phys. 43, 871 (2005)), there appears to be a difference in the exponents of a factor of two. The source code seems to indicate that the factor of 2 is used in the computation. Just wanted to check before fitting the wrong functional form!

mkanski · November 18, 2022, 12:16pm

This appears to be a typo in the documentation, the factor of 2 is in the source code:

github.com

lammps/lammps/blob/94cc3f65900b3b64f80b360c36a22eabce4474ee/src/EXTRA-MOLECULE/bond_gaussian.cpp#L95


      
          delz = x[i1][2] - x[i2][2];
          
          
rsq = delx * delx + dely * dely + delz * delz;
          r = sqrt(rsq);
          
          
sum_g_i = 0.0;
          sum_numerator = 0.0;
          for (int i = 0; i < nterms[type]; i++) {
            dr = r - r0[type][i];
            prefactor = (alpha[type][i] / (width[type][i] * sqrt(MY_PI2)));
            exponent = -2 * dr * dr / (width[type][i] * width[type][i]);
            g_i = prefactor * exp(exponent);
            sum_g_i += g_i;
            sum_numerator += g_i * dr / (width[type][i] * width[type][i]);
          }
          
          
// force & energy
          if (sum_g_i < SMALL) sum_g_i = SMALL;
          
          
if (r > 0.0)
            fbond = -4.0 * (force->boltz * bond_temperature[type]) * (sum_numerator / sum_g_i) / r;

Same for angle:

github.com

lammps/lammps/blob/94cc3f65900b3b64f80b360c36a22eabce4474ee/src/EXTRA-MOLECULE/angle_gaussian.cpp#L126


      
          s = 1.0 / s;
          
          
// force & energy
          double theta = acos(c);
          
          
sum_g_i = 0.0;
          sum_numerator = 0.0;
          for (int i = 0; i < nterms[type]; i++) {
            dtheta = theta - theta0[type][i];
            prefactor = (alpha[type][i] / (width[type][i] * sqrt(MY_PI2)));
            exponent = -2 * dtheta * dtheta / (width[type][i] * width[type][i]);
            g_i = prefactor * exp(exponent);
            sum_g_i += g_i;
            sum_numerator += g_i * dtheta / (width[type][i] * width[type][i]);
          }
          
          
if (sum_g_i < SMALL) sum_g_i = SMALL;
          if (eflag) eangle = -(force->boltz * angle_temperature[type]) * log(sum_g_i);
          
          
// I should check about the sign of this expression
          a = -4.0 * (force->boltz * angle_temperature[type]) * (sum_numerator / sum_g_i) * s;

akohlmey · November 18, 2022, 6:25pm

I have checked the output from bond_write for bond style gaussian against the expression in the docs and can confirm that the factor 2 is missing (and a right parenthesis too many).
This will be updated with the next LAMMPS patch release (probably next month).

Qijun_1 · November 22, 2022, 2:27am

Hi,I have meet the same problem these days, but I don’t know whether my
problem is because of this missing factor 2, I compared the gaussian type (n=1) and transformed it to harmonic type ‘K(l-l0)^2’, ignored energy when L=L0, but bond length in gaussian simulation have smoe much longer bond length as 90A rather than 4.8A in harmonic simulation. Could you
please help me?

else I also see that original paper, I still don’t konw why they normalized by P/r2 or P/sin(theta)

akohlmey · November 22, 2022, 3:21am

Please see my response here Bond_style Gaussian question in simulation (n=1,compared with harmonic) - #6 by akohlmey