Error for angle_style gaussian

[xih72]

Dear all,

I observed unexpected wrong results from “angle_style gaussian”, while the results from “bond_style gaussian” is correct.

I tested a bulk system of coarse-grained (CG) decane molecules. Each CG decane molecule contains 5 CG atoms. Each CG atom is either Cy type representing CH3CH2- or CE type representing -CH2CH2-.

Based on all-atom (AA) simulation, the distributions of pseudo Cy-CE & CE-CE bonds have two peaks, and the distributions of pseudo Cy-CE-CE & CE-CE-CE angles have 4 peaks (one is weak). I want to reproduce these distributions with a CG model, therefore, I decide to adopt the Gaussian Sum functional form for CG bonds and angles:

units real
atom_style full
bond_style gaussian
angle_style gaussian
dihedral_style charmm
pair_style lj/charmm/coul/long 14.0 18.0 13.2
special_bonds amber
kspace_style ewald 1e-6
read_data data.bulk

In the data.bulk file, I set two sets of gaussian term parameters for each bond type
Bond Coeffs
1 298.15 2 0.401 0.120 2.330 0.601 0.102 2.581
2 298.15 2 0.342 0.120 2.337 0.660 0.101 2.587
And the simulation of CG system correctly lead to a two-peak distribution for each CG bond type, which perfectly match the result from AA simulations.

However, the resulted distribution for each CG angle type exhibits only one same broad peak, no matter whether I set two, three, or four sets of gaussian term parameters for each CG angle type in the data.bulk file:
Angle Coeffs
1 298.15 2 0.398 11.105 154.834 0.274 8.113 173.860
2 298.15 2 0.377 11.100 155.018 0.292 8.009 173.907
OR:
1 298.15 3 0.226 10.338 126.681 0.398 11.105 154.834 0.274 8.113 173.860
2 298.15 3 0.239 9.984 126.792 0.377 11.100 155.018 0.292 8.009 173.907
OR:
1 298.15 4 0.114 10.344 112.049 0.226 10.338 126.681 0.398 11.105 154.834 0.274 8.113 173.860
2 298.15 4 0.105 10.125 112.218 0.239 9.984 126.792 0.377 11.100 155.018 0.292 8.009 173.907

The version of LAMMPS which I use is: LAMMPS (29 Aug 2024 - Update 2)
Therefore, something must be wrong with the “angle_style gaussian” in this version.
Can anyone give a suggestion on how to debug this issue? Which source file should I check?
Thanks!

[akohlmey]

@xih72

Therefore, something must be wrong with the “angle_style gaussian” in this version.

That is unlikely.

It is more likely that your choice of coefficients does not yield the potential you are expecting.

@xih72

Which source file should I check?

I don’t think that the issue is in the source code.

You can easily check what potential energy (and force) you get by using the angle_write command. It should be sufficient to check the potential energy because we check with fix numdiff that the forces and the potential energy match.

Here is a minimal example input using one of your parameter sets:

units real
atom_style angle

region box block -5 5 -5 5 -5 5
create_box 1 box angle/types 2
mass 1 1.0

angle_style gaussian
angle_coeff 1 298.15 3 0.226 10.338 126.681 0.398 11.105 154.834 0.274 8.113 173.860
angle_coeff 2 298.15 3 0.239 9.984 126.792 0.377 11.100 155.018 0.292 8.009 173.907

shell rm -f angle_table.dat
angle_write 1 100 angle_table.dat Gaussian1
angle_write 2 100 angle_table.dat Gaussian2

I am attaching 3 images, one plotting your first set of parameters, one with your second set of parameters, and a third replacing the angle type 1 parameters with the ones from the example in the documentation. The difference is quite telling.

[xih72]

Hi Axel,

Thank you very much for your reply!

I did some work based on your suggestion. Here are what I have found:

The example parameters of angle_coeff for “angle_style gaussian” in the documentation are obviously copied from the example of “bond_style gaussian”, but it does give me a hint that the wi values are too large leading to very broad peaks.

The sets of parameters given in my first post are from least-square fit with the Levenberg-Marquardt method of AA simulations, as the paper by Milano et al. (2005) described. The magnitudes of the parameter values also agree with the Table 3 in Milano’s paper. But I do notice that the documentation of “angle_style gaussian” mentions that the units of Ai & wi are radians. Therefore, I re-calculated the distributions of pseudo angles from AA simulations as a function of radians [0, pi] instead of degrees [0, 180], and the least-square fit. Now all wi parameters are reduced by ~ 57 (=180.0/pi) fold.
angle_coeff 1 298.15 4 0.1150 0.1810 112.098 0.2250 0.1790 126.691 0.4030 0.1980 154.900 0.2590 0.1300 173.754
angle_coeff 2 298.15 4 0.1060 0.1770 112.254 0.2380 0.1730 126.795 0.3830 0.1980 155.096 0.2750 0.1280 173.807

But we can notice that the values of Ai terms are still the same as before. The angle distributions of CG simulations do have separated peaks at the same peak positions of AA result (see attached picture)

Unfortunately the amplitudes (and hence the areas) of the peaks do not fit perfectly as expected. (In Milano’s paper, “no automatic optimization was needed to adjust the potential parameter Ai, wi, and lci or Thetaci. The values obtained from the simple Boltzmann inversion of Gaussian-fitted target distributions already reproduced well the peak position, width, and area.”)

The amplitudes (and hence the areas) of the peaks can be optimized by further adjusting the Ai parameters. I noticed that if Ai values are set as small as the order shown in the documentation example, the simulations are not so stable (easy to crash).

Also, I noticed that the “angle_write” command can only be executed with single CPU core. If multiple CPU cores are used, the LAMMPS (29 Aug 2024 - Update 2) would crash with the following error message:

Creating table file angle_table.dat with DATE: 2025-07-07
[qin:53552] *** An error occurred in MPI_Comm_free
[qin:53552] *** reported by process [1036386305,4613425729290371073]
[qin:53552] *** on communicator MPI_COMM_WORLD
[qin:53552] *** MPI_ERR_COMM: invalid communicator
[qin:53552] *** MPI_ERRORS_ARE_FATAL (processes in this communicator will now abort,
[qin:53552] *** and potentially your MPI job)

[akohlmey]

@xih72

Also, I noticed that the “angle_write” command can only be executed with single CPU core. If multiple CPU cores are used, the LAMMPS (29 Aug 2024 - Update 2) would crash with the following error message:

Thanks for the report. This is easily fixed. See: Collected small changes and fixes by akohlmey · Pull Request #4650 · lammps/lammps · GitHub