[gmx-users] pull=constraint gives zero forces

Thomas Schlesier schlesi at uni-mainz.de
Tue Oct 9 15:04:24 CEST 2012


For me it seems that the problem with the constraints iteration is only 
relevant if both groups are connected with contraints over the rest of 
the system.
It seems that your pull-code parameters are essentially the same like 
mine (apart from the fact that you have more groups).
In my case the system consisted of a dimer, and the pull-constraint 
acted between both monomers (which had internal constraints) and i had 
never any problems when the seperation between the two molecules was 
reasonable.

But nevertheless if the system breaks if you exchange the constraints to 
normal bonds, the is somewhere other a (probably additional) problem.
Is the system stable, when you do not use the pull-code?

Am 09.10.2012 14:26, schrieb gmx-users-request at gromacs.org:
> Hi Erik,
>
> Yes - I had a feeling I'd read that in the manual somewhere, which was
> why I tried replacing such constraints with bonds. As I wrote above,
> the LINCS warnings disappeared but the segfault remained.
>
> Alex
>
> 2012/10/9 Erik Marklund [via GROMACS]<ml-node+s5086n5001821h19 at n6.nabble.com>:
>> >Hi,
>> >
>> >Do you know there are issues with using pull=constraint on molecules that
>> >have constrained bonds? It's mentioned in the manual somewhere.
>> >
>> >Erik
>> >
>> >9 okt 2012 kl. 11.39 skrev alex.bjorling:
>> >
>>> >>Sorry - forgot to mention that before crashing, the run with all other
>>> >>constraints removed produces a single line of pullf output:
>>> >>
>>> >>0.0000  -812.401        -4002.84        482.04  1951.47 138.953 -1806.55
>>> >>-601.007        2644.79 447.018 1768.6  -214.64 -199.829        -2746.97
>>> >>1177.7  476.39  288.535 -559.274        123.08  114.493 851.86  550.558
>>> >>
>>> >>As Thomas Schlesier mentions here,
>>> >>
>>> >>http://gromacs.5086.n6.nabble.com/pull-constraint-gives-zero-forces-tp5001817.html,
>>> >>the pullf output apparently contains the forces necessary to enforce the
>>> >>constraints.
>>> >>
>>> >>/ Alex
>>> >>
>>> >>
>>> >>alex.bjorling wrote
>>>> >>>Thanks guys,
>>>> >>>
>>>> >>>Fixing the bug, recompiling and trying again results in a segfault
>>>> >>>like with version 4.0.7. I interpret this as GROMACS working fine, and
>>>> >>>suppose that there's something wrong with my input. Will continue this
>>>> >>>thread here and would really appreciate any ideas on how to proceed.
>>>> >>>
>>>> >>>Before the segfault, I get a bunch of LINCS warnings, all concerning
>>>> >>>atoms that have other constraints in the topology. Manually replacing
>>>> >>>these by stiff bonds in the itp gets rid of the LINCS warnings, but
>>>> >>>still produces an immediate segfault.
>>>> >>>
>>>> >>>The complete mdp follows. (Chris: previously posted via the web forum
>>>> >>>- it seems then you have to click the nabble link to see it).
>>>> >>>
>>>> >>>Cheers,
>>>> >>>Alex
>>>> >>>
>>>> >>>50_constr3.mdp:
>>>> >>>******************************************************************
>>>> >>>pull            = constraint
>>>> >>>pull_geometry   = position
>>>> >>>pull_dim        = Y Y Y
>>>> >>>pull_start      = yes
>>>> >>>pull_group0     = BB
>>>> >>>pull_nstxout    = 1000
>>>> >>>pull_nstfout    = 1000
>>>> >>>pull_ngroups    = 7
>>>> >>>pull_constr_tol = 1e-6
>>>> >>>
>>>> >>>pull_group1     = group1
>>>> >>>pull_init1      = 0.0 0.0 0.0
>>>> >>>pull_rate1      = 0.0 0.0 0.0
>>>> >>>
>>>> >>>pull_group2     = group2
>>>> >>>pull_init2      = 0.0 0.0 0.0
>>>> >>>pull_rate2      = 0.0 0.0 0.0
>>>> >>>
>>>> >>>pull_group3     = group3
>>>> >>>pull_init3      = 0.0 0.0 0.0
>>>> >>>pull_rate3      = 0.0 0.0 0.0
>>>> >>>
>>>> >>>pull_group4     = group4
>>>> >>>pull_init4      = 0.0 0.0 0.0
>>>> >>>pull_rate4      = 0.0 0.0 0.0
>>>> >>>
>>>> >>>pull_group5     = group5
>>>> >>>pull_init5      = 0.0 0.0 0.0
>>>> >>>pull_rate5      = 0.0 0.0 0.0
>>>> >>>
>>>> >>>pull_group6     = group6
>>>> >>>pull_init6      = 0.0 0.0 0.0
>>>> >>>pull_rate6      = 0.0 0.0 0.0
>>>> >>>
>>>> >>>pull_group7     = group7
>>>> >>>pull_init7      = 0.0 0.0 0.0
>>>> >>>pull_rate7      = 0.0 0.0 0.0
>>>> >>>
>>>> >>>;
>>>> >>>; STANDARD MD INPUT OPTIONS FOR MARTINI 2.0 or 2.1
>>>> >>>;
>>>> >>>
>>>> >>>; TIMESTEP IN MARTINI
>>>> >>>; Most simulations are numerically stable
>>>> >>>; with dt=40 fs, some (especially rings) require 20-30 fs.
>>>> >>>; Note that time steps of 40 fs and larger may create local heating or
>>>> >>>; cooling in your system. Although the use of a heat bath will globally
>>>> >>>; remove this effect, it is advised to check consistency of
>>>> >>>; your results for somewhat smaller time steps in the range 20-30 fs.
>>>> >>>; Time steps exceeding 40 fs should not be used; time steps smaller
>>>> >>>; than 20 fs are also not required.
>>>> >>>
>>>> >>>;define = -DPOSRES
>>>> >>>integrator               = md
>>>> >>>tinit                    = 0.0
>>>> >>>dt                       = 0.02
>>>> >>>nsteps                   = 2500000 ; 50 ns
>>>> >>>nstcomm                  = 1
>>>> >>>comm-grps =
>>>> >>>
>>>> >>>nstxout                  = 0
>>>> >>>nstvout                  = 0
>>>> >>>nstfout                  = 0
>>>> >>>nstlog                   = 1000
>>>> >>>nstenergy                = 100
>>>> >>>nstxtcout                = 1000
>>>> >>>xtc_precision            = 100
>>>> >>>xtc-grps                 =
>>>> >>>energygrps               = Protein
>>>> >>>
>>>> >>>; NEIGHBOURLIST and MARTINI
>>>> >>>; Due to the use of shifted potentials, the noise generated
>>>> >>>; from particles leaving/entering the neighbour list is not so large,
>>>> >>>; even when large time steps are being used. In practice, once every
>>>> >>>; ten steps works fine with a neighborlist cutoff that is equal to the
>>>> >>>; non-bonded cutoff (1.2 nm). However, to improve energy conservation
>>>> >>>; or to avoid local heating/cooling, you may increase the update
>>>> >>>frequency
>>>> >>>; and/or enlarge the neighbourlist cut-off (to 1.4 nm). The latter option
>>>> >>>; is computationally less expensive and leads to improved energy
>>>> >>>conservation
>>>> >>>
>>>> >>>nstlist                  = 10
>>>> >>>ns_type                  = grid
>>>> >>>pbc                      = xyz
>>>> >>>rlist                    = 1.4
>>>> >>>
>>>> >>>; MARTINI and NONBONDED
>>>> >>>; Standard cut-off schemes are used for the non-bonded interactions
>>>> >>>; in the Martini model: LJ interactions are shifted to zero in the
>>>> >>>; range 0.9-1.2 nm, and electrostatic interactions in the range 0.0-1.2
>>>> >>>nm.
>>>> >>>; The treatment of the non-bonded cut-offs is considered to be part of
>>>> >>>; the force field parameterization, so we recommend not to touch these
>>>> >>>; values as they will alter the overall balance of the force field.
>>>> >>>; In principle you can include long range electrostatics through the use
>>>> >>>; of PME, which could be more realistic in certain applications
>>>> >>>; Please realize that electrostatic interactions in the Martini model are
>>>> >>>; not considered to be very accurate to begin with, especially as the
>>>> >>>; screening in the system is set to be uniform across the system with
>>>> >>>; a screening constant of 15. When using PME, please make sure your
>>>> >>>; system properties are still reasonable.
>>>> >>>
>>>> >>>coulombtype              = Shift
>>>> >>>rcoulomb_switch          = 0.0
>>>> >>>rcoulomb                 = 1.2
>>>> >>>epsilon_r                = 15
>>>> >>>vdw_type                 = Shift
>>>> >>>rvdw_switch              = 0.9
>>>> >>>rvdw                     = 1.2
>>>> >>>DispCorr                 = No
>>>> >>>
>>>> >>>; MARTINI and TEMPRATURE/PRESSURE
>>>> >>>; normal temperature and pressure coupling schemes can be used.
>>>> >>>; It is recommended to couple individual groups in your system
>>>> >>>separately.
>>>> >>>; Good temperature control can be achieved with the Berendsen thermostat,
>>>> >>>; using a coupling constant of the order of ????? = 1 ps. Even better
>>>> >>>; temperature control can be achieved by reducing the temperature
>>>> >>>coupling
>>>> >>>; constant to 0.1 ps, although with such tight coupling (????? approaching
>>>> >>>; the time step) one can no longer speak of a weak-coupling scheme.
>>>> >>>; We therefore recommend a coupling time constant of at least 0.5 ps.
>>>> >>>;
>>>> >>>; Similarly, pressure can be controlled with the Berendsen barostat,
>>>> >>>; with a coupling constant in the range 1-5 ps and typical
>>>> >>>compressibility
>>>> >>>; in the order of 10-4 - 10-5 bar-1. Note that, in order to estimate
>>>> >>>; compressibilities from CG simulations, you should use Parrinello-Rahman
>>>> >>>; type coupling.
>>>> >>>
>>>> >>>tcoupl                   = Berendsen
>>>> >>>tc-grps                  = system
>>>> >>>tau_t                    = 1.0
>>>> >>>ref_t                    = 320
>>>> >>>Pcoupl                   = berendsen
>>>> >>>Pcoupltype               = isotropic
>>>> >>>tau_p                    = 2.0
>>>> >>>compressibility          = 3e-4
>>>> >>>ref_p                    = 1.0
>>>> >>>
>>>> >>>gen_vel                  = no
>>>> >>>gen_temp                 = 320
>>>> >>>gen_seed                 = 473529
>>>> >>>
>>>> >>>; MARTINI and CONSTRAINTS
>>>> >>>; for ring systems constraints are defined
>>>> >>>; which are best handled using Lincs.
>>>> >>>; Note, during energy minimization the constrainst should be
>>>> >>>; replaced by stiff bonds.
>>>> >>>
>>>> >>>constraints              = none
>>>> >>>constraint_algorithm     = Lincs
>>>> >>>unconstrained_start      = no
>>>> >>>lincs_order              = 4
>>>> >>>lincs_warnangle          = 30
>>>> >>>
>>>> >>>******************************************************************
>>>> >>>
>>>> >>>
>>>> >>>
>>>> >>>2012/10/9 Jaakko Uusitalo [via GROMACS]
>>>> >>>&lt;
>>> >>
>>>> >>>ml-node+s5086n5001810h78 at .nabble
>>> >>
>>>> >>>&gt;:
>>>>> >>>>
>>>>> >>>>constraint pulling has a bug in 4.5.5, see:
>>>>> >>>>http://redmine.gromacs.org/issues/825. I'm guessing that's causing your
>>>>> >>>>problems. Fixing it is very easy (see the link) or you can also use an
>>>>> >>>>earlier version like 4.5.3 that works.
>>>>> >>>>
>>>>> >>>>On 9.10.12 2:26 , Christopher Neale wrote:
>>>>> >>>>
>>>>>> >>>>>Please post your entire .mdp file and a snip of the output in your
>>>>>> >>>>>pullf
>>>>>> >>>>>and pullc files.
>>>>>> >>>>>(Your initial post on this topic was also missing these, although the
>>>>>> >>>>>text
>>>>>> >>>>>reads as if you intended to include them).
>>>>>> >>>>>
>>>>>> >>>>>Chris.
>>>>>> >>>>>
>>>>>> >>>>>-- original message --
>>>>>> >>>>>
>>>>>> >>>>>Following up on this post. I've tried the same runs using version
>>>>>> >>>>>4.0.7,
>>>>>> >>>>>which gave immediate segmentation faults. Not sure if this is a clue or
>>>>>> >>>>>a
>>>>>> >>>>>trivial consequence of switching versions, but there it is.
>>>>>> >>>>>
>>>>>> >>>>>Any other ideas why the pullf output just contains zeros?
>>>>>> >>>>>
>>>>>> >>>>>Cheers,
>>>>>> >>>>>Alex
>>> >>




More information about the gromacs.org_gmx-users mailing list