Inquiry Regarding the Implementation of Angle Function Type 4 of GROMACS in setup_gmx code

[Dropletsimuli]

Hi, I am attempting to set up a GROMACS simulation with SPICA using setup_gmx, but I am having some difficulty in understanding the specifics of how function type 4 (cross bond-angle potential) is implemented. According to the GROMACS documentation(https://manual.gromacs.org/current/reference-manual/topologies/topology-file-formats.html#), function type 4 is used to represent a bond-angle potential with cross terms(https://manual.gromacs.org/current/reference-manual/functions/bonded-interactions.html#bond-angle-cross-term), described by the formula:

However, in the code, I noticed that the parameters defined for function type 4 (tmpcalca, tmpcalcb, sig, eps) do not seem to directly correspond to the parameters in the above formula(ri, rj, rk, and 𝑘𝑟𝜃) . In the setup_gmx code:

tmpcalca = database.ange[ang_params[uniq_angs-1]] tmpcalcb = database.fang[ang_params[uniq_angs-1]] 4.184 2.0 eps = database.eps[vdwtmp]*4.184 sig = database.sig[vdwtmp]/10.0

Given this context, I would like to ask the following questions:

1. In the implementation of function type 4, how do the parameters tmpcalca, tmpcalcb, sig, and eps correspond to the parameters defined for the cross bond-angle potential in GROMACS?
2. If these parameters do not correspond to the standard function type 4, what role do they play in the simulation?
3. Is there any specific adjustment I need to make to ensure the correct implementation of function type 4?

Your guidance would be greatly helpful in advancing my simulation work. Any insights or guidance you can provide would be greatly appreciated.

Best regards, Drep

[yskmiyazaki]

@Dropletsimuli Thank you for trying to use SPICA with GROMACS.

We need to use an angle potential that is not implemented in GROMACS by default. So, for MD simulations with SPICA, we have to apply a patch file to replace the cross bond-angle potential with the SPICA angle potential. See https://github.com/SPICA-group/gromacs-spica.

After correctly applying the patch, you can use the following angle potential used in SPICA as "type 4" instead of the cross bond-angle potential. $U{angle} = U{1-2-3} + U{1-3}$ $U{1-2-3} = k{\theta} (\theta{ijk} - \theta{0})$ $U{1-3} = U{LJ}(r{ij}; \varepsilon{ij}, \sigma{ij}) - U{LJ}(\sigma{ij})\ for\ r{ij} < \sigma{ij}$ This is a combination of a harmonic potential term and a replusive term of an LJ potential. See https://www.spica-ff.org/forcefield.html for more details.

As you mentioned, the cross bond-angle potential has four parameters $r{1e}, r{2e}, r{3e},\ and\ k{r\theta}$. After applying the patch, the first - fourth arguments correspond to the following parameters in the SPICA angle potential: i. $\theta0$ : equilibrium angle ii. $k{\theta}$ : force constant
iii. $\sigma$ : LJ parameter iv. $\varepsilon$ : LJ parameter

With these things in mind, the answers to your questions are as follows:

1. The parameters found in setup_gmx tmpcalca tmpcalcb, eps, sig correspond to $\theta0$, $k{\theta}$, $\varepsilon$, and $\sigma$ in the SPICA angle potential. So, they do not correspond to the cross bond-angle term parameters.
2. After applying the patch file, they have meanings as parameters defined in the SPICA potential as mentioned above.
3. Yes, there is a patch file to replace the "type 4" angle function to the SPICA potential. Please have a look at https://github.com/SPICA-group/gromacs-spica.

Should you still have any queries after reading my answers, please feel free to ask me.

[Dropletsimuli]

@Dropletsimuli Thank you for trying to use SPICA with GROMACS.

We need to use an angle potential that is not implemented in GROMACS by default. So, for MD simulations with SPICA, we have to apply a patch file to replace the cross bond-angle potential with the SPICA angle potential. See https://github.com/SPICA-group/gromacs-spica.

After correctly applying the patch, you can use the following angle potential used in SPICA as "type 4" instead of the cross bond-angle potential. Uangle=U1−2−3+U1−3 U1−2−3=kθ(θijk−θ0) U1−3=ULJ(rij;εij,σij)−ULJ(σij) for rij<σij This is a combination of a harmonic potential term and a replusive term of an LJ potential. See https://www.spica-ff.org/forcefield.html for more details.

As you mentioned, the cross bond-angle potential has four parameters r1e,r2e,r3e, and krθ. After applying the patch, the first - fourth arguments correspond to the following parameters in the SPICA angle potential: i. θ0 : equilibrium angle ii. kθ : force constant iii. σ : LJ parameter iv. ε : LJ parameter

With these things in mind, the answers to your questions are as follows:

1. The parameters found in setup_gmx tmpcalca tmpcalcb, eps, sig correspond to θ0, kθ, ε, and σ in the SPICA angle potential. So, they do not correspond to the cross bond-angle term parameters.
2. After applying the patch file, they have meanings as parameters defined in the SPICA potential as mentioned above.
3. Yes, there is a patch file to replace the "type 4" angle function to the SPICA potential. Please have a look at https://github.com/SPICA-group/gromacs-spica.

Should you still have any queries after reading my answers, please feel free to ask me.

I got it. This will be a boost to my progress in our work. I deeply appreciate your assistance and expertise. Looking forward to possibly collaborating again in the future!

Thanks again! @yskmiyazaki