Transformation pull co-ordinate dihedral maths is odd

User discussions
,
1

GROMACS version: 2022.3
GROMACS modification: No

I have these transformations pull co-ordinate for dihedral in the mdp file for 4-Ala as a test case.

pull-coord5-geometry      = transformation
pull-coord5-groups         = []
pull-coord5-type                = external-potential
pull-coord5-potential-provider  = AWH
pull-coord5-expression = 0.5*x1+0.5*x3
;
pull-coord6-geometry            = transformation
pull-coord6-groups              = []
pull-coord6-type                = external-potential
pull-coord6-potential-provider  = AWH
pull-coord6-expression          = 0.5*x2+0.5*x4
;
awh                      = yes
awh-potential            = convolved
awh-share-multisim       = yes
awh-nbias                = 1
awh-nstout               = 50000
awh1-ndim                    = 2
awh1-equilibrate-histogram   = yes
awh1-target                  = constant
awh1-share-group             = 1
awh1-dim1-coord-index        = 5
awh1-dim2-coord-index        = 6

awh1-dim1-start          = -179
awh1-dim1-end            = 179
awh1-dim1-diffusion      = 2e-4
awh1-dim1-force-constant  = 10

awh1-dim2-start          = -179
awh1-dim2-end            = 179
awh1-dim2-diffusion      = 2e-4
awh1-dim2-force-constant  = 10

and these are my pullx file output.

each column - the first column is time, the first is phi1, the second is psi1, third is phi2, fourth is psi2 and the last two is 0.5x1+0.5x3, 0.5x2+0.5x4.

Without being good at maths, I’m sure that (-37.701±52.4284)/2 is not equal to -0.786528. Not in radian - and maybe not as a sin or cos function.

Is this a bug, or am I missing a feature?

0.0000  -37.701 -37.701 -52.4284        26.4062 -0.786528       -0.0985655
0.1000  -30.6594        -30.6594        -41.1388        21.624  -0.626558       -0.0788487
0.2000  -33.2508        -33.2508        -39.729 21.0787 -0.636869       -0.106222
0.3000  -30.0481        -30.0481        -36.201 19.2426 -0.578132       -0.0942957
0.4000  -24.32  -24.32  -41.1157        16.7663 -0.571034       -0.0659185
0.5000  -33.3076        -33.3076        -32.1182        22.0187 -0.570947       -0.0985144
0.6000  -29.0465        -29.0465        -48.9284        27.1546 -0.68046        -0.0165097
0.7000  -24.6533        -24.6533        -46.6624        20.2891 -0.622348       -0.0380849
0.8000  -14.0262        -14.0262        -40.0957        8.97931 -0.472303       -0.0440426
0.9000  -30.801 -30.801 -44.182 22.107  -0.654351       -0.0758701
1.0000  -32.0729        -32.0729        -42.2959        29.7339 -0.648991       -0.0204121
1.1000  -13.0682        -13.0682        -48.9984        32.2253 -0.541634       0.167177
1.2000  -21.3409        -21.3409        -41.3596        32.7632 -0.547165       0.0996788

Best wishes

Will

2

After some mathematical decryption in the morning, I have now found out that the angles in degree is converted to radian for the transformation pull co-ordinate calculation and NOT degrees. The output of each term - however, is in degrees, but the output from the maths - is in radian.

The answer is

-37.701 → -0.658
-52.4284 → -0.9150

(-0.658 + -0.9150)/2 = -0.7865

which then matches the 5th column.

This should be noted in the documentation for clarity in the future.

3

However, during the simulation, x1 only sampled from -60 to +60 degree. The bound of the sum should let it sample from -180 to +180 given the bound is now set to -3.1 to 3.1.

Any idea?

4

In addition to above, I have conducted several test cases as follow:

  1. Just doing AWH on one dihedral (phi) on the Alax4 with the bound at -180 to -100. It still only samples between -30 and -60.

This suggested that the problem is not only at the transformation pull co-ordinate, but at the dihedral implementation itself.

pull-coord1-geometry        = dihedral
pull-coord1-groups          = 3 1 1 2 2 4
pull-coord1-type                = external-potential
pull-coord1-potential-provider  = AWH


pull-coord2-geometry        = dihedral
pull-coord2-groups          = 4 2 2 1 1 3

pull-coord3-geometry        = dihedral
pull-coord3-groups          = 6 4 4 5 5 7

pull-coord4-geometry        = dihedral
pull-coord4-groups          = 7 5 5 6 6 8

;pull-coord5-geometry      = transformation
;pull-coord5-groups         = []
;pull-coord5-type                = external-potential
;pull-coord5-potential-provider  = AWH
;pull-coord5-expression = 0.5*x1+0.5*x3
;
;pull-coord6-geometry            = transformation
;pull-coord6-groups              = []
;pull-coord6-type                = external-potential
;pull-coord6-potential-provider  = AWH
;pull-coord6-expression          = 0.5*x2+0.5*x4
;
awh                      = yes
awh-potential            = convolved
awh-share-multisim       = yes
awh-nbias                = 1
awh-nstout               = 50000
awh1-ndim                    = 1
awh1-equilibrate-histogram   = yes
awh1-target                  = constant
awh1-share-group             = 1
awh1-dim1-coord-index        = 1
;awh1-dim2-coord-index        = 6

awh1-dim1-start          = -180
awh1-dim1-end            = -100
awh1-dim1-diffusion      = 0.0001
awh1-dim1-force-constant  = 10

;awh1-dim2-start          = -3.1
;awh1-dim2-end            = 3.1
;awh1-dim2-diffusion      = 20
5

I currently see no reason, either, why the sampling should be limited and not cover the full range from -180 to +180 degrees or the range you specified. However:

  1. At which dihedral angle did you start your simulations?
  2. How long are these simulations? Did AWH have a chance to build up the bias necessary?
6

The initial simulation starts at -35 deg on the phi angle.

I had cases where I ran for 200 ns - and those region hasn’t been explored or covered. I would suspect it to be quite fast with alanine quadpeptide though?

7

Sure, it’s a small system, but as far as I remember, it was shown that one dihedral alone is a bad choice for the RC such that you might just be unlucky and still be stuck in the initial free-energy well. I’d say it’s not a clear indicator of a bug because poor sampling is still a possible explanation.

If you’re really worried about a wrong dihedral implementation, test alanine dipeptide in vacuum, use both backbone dihedrals as RC and see if you can sample the entire range.

8

The force constant for a dihedral coordinate is in kJ/mol/rad^2, as noted in the use guide. Your value of 10 seems rather small to me. Likely it’s too small to force the coordinate to sample the whole range. You can see if this is the case by comparing the coordinate distribution and the reference distribution in the output of gmx awh -more

I pushed up a note to the user guide that angular coordinates use units of radians in a transformation coordinate.