[gmx-users] constraining multiple types of bonds

[tarak karmakar]

Dear All,


In my system I need to fix three types of bonds

1) Metal-Ligand distance at a particular value given in PDB ( not covalent)
2) I need to fix some of the bond lengths (covalent) for the substrate molecule.
3) Lastly the covalent H-bonds ( C-H, N-H, O-H etc.)


My input .mdp file is given below

; 7.3.3 Run Control
integrator              = md-vv                 ; md integrator
tinit                   = 0                         ; [ps] starting time for run
dt                      = 0.001                 ; [ps] time step for integration
nsteps                  = 5000000               ; maximum number of
steps to integrate, 0.001 * 20,00,000 = 2 ns
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                 = 5000                 ; [steps] freq to write
coordinates to trajectory
nstvout                 = 5000                 ; [steps] freq to write
velocities to trajectory
nstfout                 = 5000                 ; [steps] freq to write
forces to trajectory
nstlog                  = 1000                 ; [steps] freq to write
energies to log file
nstenergy               = 1000                 ; [steps] freq to write
energies to energy file
nstxtcout               = 1000                 ; [steps] freq to write
coordinates to xtc trajectory
xtc_precision           = 1000                 ; [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                   = 1.2                   ; [nm] cut-off
distance for the short-range neighbor list

nsttcouple              = 1
nstpcouple              = 1

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

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

; 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                  = Nose-Hoover                   ; Nose-Hoover
temperature coupling
tc_grps                 = Protein    Non-Protein        ; groups to
couple seperately to temperature bath
tau_t                   = 1.0        1.0                ; [ps] time
constant for coupling
ref_t                   = 300        300                ; [K]
reference temperature for coupling

; 7.3.15 Pressure Coupling
pcoupl                  = MTTK                          ; pressure
coupling where box vectors are variable
pcoupltype              = isotropic                     ; pressure
coupling in x-y-z directions
tau_p                   = 1.0                           ; [ps] time
constant for coupling
compressibility         = 4.5e-5                        ; [bar^-1]
compressibility
ref_p                   = 1.0                           ; [bar]
reference pressure for coupling

; 7.3.17 Velocity Generation
gen_vel                 = no                            ; velocity
generation turned off

; 7.3.18 Bonds
constraints             = h-bonds
constraint_algorithm    = SHAKE                         ; SHAKE
Constraint Solver
shake_tol               = 1.0e-5




So I'm bit confused how to implement constraints algorithm for these
type of problem. If I do use the above set up then it is showing
following error

Program mdrun, VERSION 4.5.5
Source code file: invblock.c, line: 79

Fatal error:
Double entries in block structure. Item 5247 is in blocks 1371 and 1370
 Cannot make an unambiguous inverse block.


Thanks
-- 
Tarak