Basically, in classical MD none of bonds can be broken. 

The "broken bonds effect" you describe may come from periodic conditions - your longer, simple system simulation supports this thesis. The error message may also shed some light on reason of your problems.

Hi Rafal,

thanks for you answer. Here is the last part of the ERR file, which also includes the termination message. The problem seems to be related to LINCS i.e. bond breaking.The molecules 2185 and 2186 corresponds to the C-alpha and H-Alpha atoms of the residues which is closest to the Linker and I parametrised this part of the peptide.

---------------------------------------------------------------------------------------------------------------------

Step 17426480, time 34853 (ps)  LINCS WARNING

relative constraint deviation after LINCS:

rms 0.000003, max 0.000032 (between atoms 10546 and 10548)

bonds that rotated more than 30 degrees:

 atom 1 atom 2  angle  previous, current, constraint length

   2185   2186   32.3    0.1080   0.1080      0.1080

Step 17426481, time 34853 (ps)  LINCS WARNING

relative constraint deviation after LINCS:

rms 0.068524, max 9.419783 (between atoms 2185 and 2186)

bonds that rotated more than 30 degrees:

 atom 1 atom 2  angle  previous, current, constraint length

   2185   2186   90.0    0.1080   1.1253      0.1080

Wrote pdb files with previous and current coordinates

WARNING: Listed nonbonded interaction between particles 2182 and 2185

at distance 2.546 which is larger than the table limit 2.256 nm.

This is likely either a 1,4 interaction, or a listed interaction inside

a smaller molecule you are decoupling during a free energy calculation.

Since interactions at distances beyond the table cannot be computed,

they are skipped until they are inside the table limit again. You will

only see this message once, even if it occurs for several interactions.

IMPORTANT: This should not happen in a stable simulation, so there is

probably something wrong with your system. Only change the table-extension

distance in the mdp file if you are really sure that is the reason.

Step 17426482, time 34853 (ps)  LINCS WARNING

relative constraint deviation after LINCS:

rms 63.562419, max 8737.694336 (between atoms 2185 and 2186)

bonds that rotated more than 30 degrees:

 atom 1 atom 2  angle  previous, current, constraint length

   2183   2184   44.8    0.0997   0.0997      0.0997

   2185   2186   80.6    1.1254 943.7790      0.1080

Wrote pdb files with previous and current coordinates

-------------------------------------------------------

Program gmx mdrun, VERSION 5.1.2

Source code file: /cluster/src/gromacs/5.1.2/gnu_gpu/src/gromacs/mdlib/constr.cpp, line: 555

Fatal error:

step 17426482: Water molecule starting at atom 33148 can not be settled.

Check for bad contacts and/or reduce the timestep if appropriate.

For more information and tips for troubleshooting, please check the GROMACS

website at http://www.gromacs.org/Documentation/Errors

First, reduce the time step and run nvt simulation. Furthermore, use output coordinate file for the npt simulation with a 0.002 time step. 

A few things to consider:

1. What were your target data for the parametrization?

2. Did you check water interactions with the new linker/residue? Are any waters getting too close to a given atom/group that lead to distortion?

3. What does the in vacuo conformational sampling look like, particularly around the linker? Is it sampling configurations that would lead to clashes in solution?

The covalent bonds may be not broken, but only extended longer  than the "cutoff of bonds" in the  Visual software (e.g. VMD or Pymol). Dose the configuration changed heavily in vaccum ? In general, the higher electric force in vacuum made the structure fluctuated moderately than in water and thus the simulation in vacuum can be lasted  to  300-400 ns. The simulation crashed in water may be result from the system is not relaxed enough. You can try follow the suggestions proposed by Jamoliddin Razzokov.

@Justin,

1. I parametrise via server SWISSPARAM and cgenff, they agree well atleat for the bonded params, so I just took those from there. For partial charges, I derived those from HF/6-631G* QM calculations (consistent with CHARMMFF). 

2. I checked the water-linker and water-protein interations, I could not see any distorsion as such. None of the water comes within 1 angstrom distance from these two residues. There are Cyclohexane and Succiniamide rings in the linker, they also don't distort throughout the simulation.

3. In vacuum, the Lipid-Linker-Peptide molecule folds and forms a bulk like structure which, I think, is expected. But I don't see any major clashes or anything for the linker, or for the modified peptide residue. I would like to know more about this from you, how would you analyse better for the steric clashes in the bilayer-peptide assembly?

I was wondering if the peptide structure got distorted because of the PBC for the lipid bilayer..? I am running 2 different simulations (one NVT and other one isotropic NPT)

---------------------------------------------------------------------------------------------------------------------

@Yang,

I can understand if the bonds are longer in VMD (because of the PBC), but still the simulation should not crash in Gromacs, because Gromacs takes care of the PBC, right?

I think the initial relaxation of the system maybe a problem, I am already running a simulation with what Jamaloiddin said, let me get back with the reults.

For CHARMM, you should be doing an MP2/6-31G* optimization (MP2/6-31+G* if it is charged) of the linker model(s) and deriving charges to match its dipole moment in the optimized geometry (overestimate by 20% or so to account for implicit condensed-phase polarization). Refine the charges using HF/6-31G* water interaction energies around each atom; the CGenFF paper has details.

Make sure you're using actual atom types from the protein and lipid parameter sets and not CGenFF. Mixing and matching all of that into a single molecule is insufficiently accurate.

Run the new molecule in water (not a bilayer) and compare the conformational ensembles between the vacuum and solvated states. The in vacuo simulation just tells you there is nothing catastrophically wrong with the topology. It does not tell you that the parameters are right.