[gmx-developers] FEP and dihedraltypes, ignoring B-topology parameters

[Berk Hess]

Oliver Lange wrote:

> Hi all,
>
> I think i found a small bug, in the toppush.c of gromacs-3.3.1. 
> (release version from website)
> I am not sure though, it might be deliberately programmed in the 
> current way, which just wasn't doing what i wanted.
>
> I wanted to include a B-topology of dihedraltypes to use the topology 
> with the FEP-code.
> To this end I edited in the section [dihedraltypes] (after writing out 
> the full topol.top with the option -pp of grompp)
> such that lines in the [dihedraltype]-section contained 12 parameters:
> a1 a2 a3 a4 type c1A-c6A c1B-c6B
> I checked that all parameters are read correctly by the routine 
> push_dihedraltype().
> However, when it comes to creating the actual dihedrals push_bond() 
> calls the function default_params().
> The latter function takes special care for dihedrals with an 
> if-statment like
> if( ftype==F_PDIHS || ftype == F_RBDIHS || ftype == F_IDIHS) {--
> when it is checking if the dihedral is in the dihedral-type list.
> However, when it comes to copying the parameters it uses the for-loop
> for (j=0; (j < nrfpB); j++) {
>   p->c[nrfp+j] = bFound ? pi->c[j] : 0.0;
> (where nrfp = 6 // nrfp = interaction_function[ftype].nrfpA;)
> This loop assigns the 6 B-topology parameters (7-12) in the created 
> dihedral (p->c) with the A-topology parameters (1-6) from the 
> dihedraltype.
> I changed that line into
> if( ftype==F_PDIHS || ftype == F_RBDIHS || ftype == F_IDIHS)
>          p->c[nrfp+j] = bFound ? pi->c[nrfp+j] : 0.0;
>       else
>          p->c[nrfp+j] = bFound ? pi->c[j] : 0.0;
> and now i get the B-topology dihedrals also carried over into the 
> topol.tpr. (as checked with gmxcheck -s1 -s2)
>
> As this is what i wanted, it fixed the problem for me.
> If this is not what the GROMACS crew wants, please ignore my email.


This would not only hold for dihedrals, but for all bonded interactions.

I guess we only thought of supplying B state parameters using different atom
types for the B-state, not via supplying bonded types with different B-state
parameters.
You can not just change the line in the call as you did, as now for 
dihedrals
with different atom types in the A and B state you would get for the B-state
the B-state parameters of the dihedraltype for the B-state atom types.

The most clear solution would be to not allow for B-state parameters
in the bonded type sections.

An other option would be to use the B-state bonded type parameters
only when the A and B state atom types are the same.
But I would say things get a little confusing here.

Berk.