## Bug #2606

### Free Energy Calculation -- Function type Fourier Dih. Not implemented in ip_pert

**Description**

I believe it is a known issue, in free energy calculation for the dihedral type **5**, it will pop up some error information complaining that Fourier type dihedral is not implemented, like;

*Fatal Error: Function type Fourier Dih. Not implemented in ip_pert*

This error is mainly related with two files,

*src/gromacs/topology/topsort.cpp (line 139)src/topology/ifunc.cpp (line 111)*

In the latter file, it has defined all the type of functions and parameters that Gromacs would use, in a way enumerating those constant-string-variables.

For this definition, in topology file, under `[dihedrals]`

directive using the calculation function type **5**, after executed **grompp** command, it by default will use Fourier type function to make binary **tpr** file.

```
gmx grompp –f *mdp –c *gro –p *top –o *tpr
```

If we use Gromacs built-in command **dump** this file, we will find some lines in key word “FOURDIHS” showing the dihedral type will be used.

```
gmx dump –s *tpr > *txt [any free plain text file]
```

However, if we use dihedral function **3**, which means **Ryckaert-Bellemans** dihedral type will be used for the calculation, the difference of it with the **Fourier Dihedral** is only the conversion difference. In Gromacs manual **table 5.5**, using the cross-link we can easily find this transformation.

If we define four coefficients of **Fourier Dihedral**;

```
V = 1/2*[ F1*(1+cos(phi)) + F2*(1-cos(2*phi)) + F3*(1+cos(3*phi)) + F4*(1-cos(4*phi)) ]
```

Because of the trigonometric function, `cos(2*phi) = 2*cos^2(phi) – 1`

, then the six **Ryckaert-Bellemans** will be calculated;

```
C0 = F2 + 1/2*(F1 + F3)
C1 = 1/2*(-F1 + 3*F3)
C2 = -F2 + 4*F4
C3 = -2*F3
C4 = -4*F4
C5 = 0
```

If we follow the top procedure using a new topology file to make **tpr** file and **dump** it, using the key-word “RBDIHS”, we will know the dihedral calculation has been changed to **Ryckaert-Bellemans** for dihedral type **3**. However, this **tpr** file can be used for the Free Energy calculation without any problems.

**Therefore, it comes a question, for Free Energy calculation, does it have calculation difference between Fourier Dihedral and Ryckaert-Bellemans for the same data but only at different mathematic conversion equations?**

As far as I know, I don’t think they have any difference at the perturbation theory, since for bonded information calculation it is merely a factor of difference scale values. Besides, in code area, it seems the return parameter `bPert`

, for the function `static gmx_bool ip_pert(int ftype, const t_iparams *ip)`

in module `topsort.cpp`

, is only a **bool** value.

If not luckily, here comes two possible solutions:

1) Update the **grampp** command output result.

In the **mdp** file, if `free energy`

is turned on, during **grompp** command process, transferring dihedral type **5** to corresponded dihedral **3** in the topology file if it has any, then using these converted values to generate binary **tpr** file. In this way, which means no matter whether **Fourier Dihedral** or **Ryckaert-Bellemans** is set in the topology file, the final calculation will always be based on the **Ryckaert-Bellemans**.

2) Add **bool** return value for function `static gmx_bool ip_pert`

(line 48) in `topsort.cpp`

source file, as well as some other codes for the Fourier Dihedral perturbations.

**NOTE**

For the attached files, the same **gro** file and **mdp** file are used for **mdrun**, the difference is the **top** file, for these two **top** files, the differences are at [dihedrals] directives, in where different function types are selected.