I have some questions when I using lammps. I hope to simulate course grain bead spring model using lammps. However bond_style gaussian confused me. I first use gaussian function to represent Bond distribution. The gaussian function can transformed to Harmonic if factor n=1 and ignored the energy in balanced position. Just as followed picture. In simulation, unit is set in ‘real’, ‘npt’ ensemble ,298K and 1atm.
I compared the follows two.
bond_style harmonic
Bond Coeffs
1 20.474 4.81 # A-A
&&&&&&&&&&&&&&
bond_style gaussian
Bond Coeffs # gaussian
#Bond_type T(K) n a1 w1 L1
1 298 1 0.71 0.17020 4.808 # A-A
But the results are much more different. In gaussian style, some of calculate bond length is much longer than the balanced length which is 90A, compared with 4.81A
In harmonic style, there is no such result, and successfully represent the length distribution.
This confused me, could you please help me solve this problem, correct the function style and tell me the reason?
If you check out this post, the documentation for bond_style gaussian is missing a factor of 2 in the exponent (i.e., it should be \exp(-2(r-r_i)^2/w_i)). Though, I’m not sure if that accounts for your huge difference.
That factor of 2 is not the (main) problem that you are running into.
The implementation of bond style gaussian (and similarly angle style gaussian) has to avoid an overflow or underflow of the term inside logarithm so that the logarithm can be correctly computed.
This, however, is done in a way that only values relatively close to r0 are allowed and for others the term inside the logarithm is capped. Once the distance is too large the force will be truncated to zero.
I am currently looking into changing the truncation to allow a larger variety of r.
