[gmx-users] Constraints not working in pull code (sometimes, sometimes not)

[chris.neale at utoronto.ca]

I agree that Justin is probably correct, although constraints should  
technically work just fine with a highly dynamic reference group. Any  
problems should show up as the system blowing up, but not in correctly  
setting the position. I think that the problems can arise, however,  
when you have multiple competing constraints. You might, for example,  
try without constraining the octane bonds -- I think that you should  
get perfect match in that case. Also, it is worth comparing in single  
and double precision.

Still, I am not sure that you actually did a proper comparison for the  
following 2 reasons:

1. you have apparently different masses for the pulled group:
> Pull group 1:    15 atoms, mass   194.194
> Pull group 1:    15 atoms, mass   146.146

2. you used very different starting depth offsets: -2.81855 is not  
very close to -1.60929

PS: I would personally prefer if you avoided jumping directly to a bug  
report unless you are sure of it. There is always at least the  
possibility that you are not using the code correctly.

Chris.

> Dear all,
>
> we are trying to simulate molecule in specified depth in octanol to
> calculate logP.
> for this purpose we have used box with slab of octanol, molecule in
> specified distance in z-axis and  this mdp:
>
> integrator      = md
> dt              = 0.002
> tinit           = 0
> nsteps          = 5125000       ; 10.250 ns
> nstcomm         = 1
> comm_mode       = Linear
> ; Output parameters
> nstxout         = 0     ; every 1 ns
> nstvout         = 0
> nstfout         = 0
> nstxtcout       = 500           ; every 1 ps
> nstenergy       = 500
> ; Bond parameters
> constraints             = all-bonds
> ; Single-range cutoff scheme
> nstlist         = 5
> ns_type         = grid
> rlist           = 1.4
> rcoulomb        = 1.4
> rvdw            = 1.4
> ; PME electrostatics parameters
> coulombtype     = PME
> fourierspacing  = 0.12
> fourier_nx      = 0
> fourier_ny      = 0
> fourier_nz      = 0
> pme_order       = 4
> ewald_rtol      = 1e-5
> optimize_fft    = yes
> ; Berendsen temperature coupling is on in two groups
> Tcoupl          = V-rescale
> tc_grps         = System
> tau_t           = 0.1
> ref_t           = 310
> ; Pressure coupling is on
> Pcoupl          = berendsen
> pcoupltype      = anisotropic   ; anisotropic pressure coupling
> tau_p           = 10
> compressibility = 0 0 4.5e-5 0 0 0  ; since octanol tries to separate from
> water and shortening xy plane
> ref_p           = 1.0 1.0 1.0 0 0 0
> ; Pull code
> pull            = constraint    ;
> pull_geometry   = distance      ; Pull along the vector connecting the two
> groups. Components can be selected with pull_dim.
> pull_dim        = N N Y         ; put potential only to axis z
> pull_start      = yes           ; define initial COM distance > 0
> pull_ngroups    = 1             ; one group to pull
> pull_group0     = OCT           ; reference group (can be empty to set it to
> 0,0,0)
> pull_group1     = DRG
> pull_rate1      = 0             ; 0.01 nm per ps = 10 nm per ns = 10 m per s
> pull_k1         = 2000          ; kJ mol^-1 nm^-2
>
> however when I list pull-x.xvg I obtain:
> @ s0 legend "0 Z"
> @ s1 legend "1 dZ"
> 0.0000  6.06871 -2.81855
> 0.0200  6.07793 -2.81861
> 0.0400  6.08787 -2.81856
> 0.0600  6.08462 -2.81855
> 0.0800  6.10028 -2.82007
>
> The situation is even more strange when almost the same mdp with DOPC
> membrane works fine as can be seen in pull-x.xvg:
> 0.0000  3.72735 -1.60929
> 0.0200  3.72727 -1.60929
> 0.0400  3.7272  -1.60929
> 0.0600  3.72716 -1.60929
> 0.0800  3.72715 -1.60929
> 0.1000  3.72717 -1.60929
>
> Differences in the working case in mdp are only this:
> compressibility = 4.5e-5 4.5e-5 4.5e-5 0 0 0
>
> ref_p           = 1.0 1.0 1.0 0 0 0
> pull_group0     = DOPC
>
> Strangerly enough - the same happens in GMX4.0.7 and in GMX4.5.3
>
> I wonder where the difference is, as both logs states that
> Will apply constraint COM pulling in geometry 'distance'
> between a reference group and 1 group
> Pull group 0:  6912 atoms, mass 100624.859
> Pull group 1:    15 atoms, mass   194.194
>
> Will apply constraint COM pulling in geometry 'distance'
> between a reference group and 1 group
> Pull group 0:  9570 atoms, mass 124631.453
> Pull group 1:    15 atoms, mass   146.146
>
> Where should I look and repair?
>
> Thank you in advance
> Karel
>