[gmx-users] constraints algorithm

[tarak karmakar]

Dear All,

     In a protein simulation I imposed one bond constraint by
incorporating the "[constraints]" block in topology file. On the other
hand in .mdp file I have imposed constraints for the covalent hydrogen
bonds by implementing the LINCS algorithm. After a short equilibration
run I see the bond length is fixed [served my purpose], but then what
algorithm it's making use of to constraint the specified bond ?
The .mdp file is given below

Thanks.

......................................................................................................................................................................................
; 7.3.2 Preprocessing
;define                  = -DPOSRES      ; defines to pass to the preprocessor

; 7.3.3 Run Control
integrator              = md                    ; md integrator
tinit                   = 0                     ; [ps] starting time for run
dt                      = 0.001                ; [ps] time step for integration
nsteps                  = 50000                ; maximum number of
steps to integrate, 0.001 * 50,000 = 500 ps
comm_mode               = Linear                ; remove center of
mass translation
nstcomm                 = 1                     ; [steps] frequency of
mass motion removal
comm_grps               = Protein Non-Protein   ; group(s) for center
of mass motion removal

; 7.3.8 Output Control
nstxout                 = 1000         ; [steps] freq to write
coordinates to trajectory
nstvout                 = 1000         ; [steps] freq to write
velocities to trajectory
nstfout                 = 1000         ; [steps] freq to write forces
to trajectory
nstlog                  = 100           ; [steps] freq to write
energies to log file
nstenergy               = 100           ; [steps] freq to write
energies to energy file
nstxtcout               = 100           ; [steps] freq to write
coordinates to xtc trajectory
xtc_precision           = 100           ; [real] precision to write
xtc trajectory
xtc_grps                = System        ; group(s) to write to xtc trajectory
energygrps              = System        ; group(s) to write to energy file

; 7.3.9 Neighbor Searching
nstlist                 = 1             ; [steps] freq to update neighbor list
ns_type                 = grid          ; method of updating neighbor list
pbc                     = xyz           ; periodic boundary conditions
in all directions
rlist                   = 0.8           ; [nm] cut-off distance for
the short-range neighbor list

; 7.3.10 Electrostatics
coulombtype             = PME           ; Particle-Mesh Ewald electrostatics
rcoulomb                = 0.8           ; [nm] distance for Coulomb cut-off

; 7.3.11 VdW
vdwtype                 = cut-off       ; twin-range cut-off with
rlist where rvdw >= rlist
rvdw                    = 0.8           ; [nm] distance for LJ cut-off
DispCorr                = EnerPres      ; apply long range dispersion
corrections

; 7.3.13 Ewald
fourierspacing          = 0.12          ; [nm] grid spacing for FFT
grid when using PME
pme_order               = 4             ; interpolation order for PME, 4 = cubic
ewald_rtol              = 1e-5          ; relative strength of
Ewald-shifted potential at rcoulomb

; 7.3.14 Temperature Coupling
tcoupl                  = Berendsen                     ; temperature coupling
tc_grps                 = Protein    Non-Protein        ; groups to
couple seperately to temperature bath
tau_t                   = 0.1        0.1                ; [ps] time
constant for coupling
ref_t                   = 300        300                ; [K]
reference temperature for coupling

; 7.3.17 Velocity Generation
gen_vel                 = yes           ; generate velocities
according to Maxwell distribution of temperature
gen_temp                = 300           ; [K] temperature for Maxwell
distribution
gen_seed                = -1            ; [integer] used to initialize
random generator for random velocities

; 7.3.18 Bonds
constraints             = h-bonds       ; constraining covalent h-bonds
constraint_algorithm    = LINCS         ; LINear Constraint Solver
continuation            = no            ; no = apply constraints to
the start configuration
lincs_order             = 4             ; highest order in the
expansion of the contraint coupling matrix
lincs_iter              = 1             ; number of iterations to
correct for rotational lengthening
lincs_warnangle         = 30            ; [degrees] maximum angle that
a bond can rotate before LINCS will complain


-- 
Tarak Karmakar
Molecular Simulation Lab.
Chemistry and Physics of Materials Unit
JNCASR