Mark is right, no two ways about it. For initial equilibration and assessing preexisting structural strains try vacuum, _much_ smaller timesteps and possibly low temperatures in vacuum, only then transfer to solvent, etc. Algorithmically, LINCS requires convergence and you already are using a pretty high LINCS order... From what I see, dt = 15 fs at 310K looks like a cowboy mode simulation in this case.
Alex On Thu, Dec 15, 2016 at 2:32 PM, Mark Abraham <mark.j.abra...@gmail.com> wrote: > Hi, > > If a simulation isn't stable with a small time step (as I think you are > saying) then moving to a larger time step is guaranteed to make that worse. > Try an even smaller time step, for a long time, and see what happens. Or > take a subset of your protein and see what happens. Or simulate in vacuo > for a while. Your topology could be unsuited to your starting structure, > e.g. some part is under a lot of tension that gets released at some point > and no finite time step can in practice deal with the velocity of the > recoil... > > Mark > > On Thu, 15 Dec 2016 23:06 Nash, Anthony <a.n...@ucl.ac.uk> wrote: > > > Hi all, > > > > I¹m hoping for some help. I¹m very sorry, this is a bit of a long one. > > > > I¹ve been struggling for almost a month trying to run a CG representation > > of our all-atom model of a collagen protein (3 polypeptide chains in a > > protein). Our original AMBER all-atom model had been successful modelling > > using MD. I went on to use the latest version of Martinize.py with the > > latest version of the MARTINI forcefield fields. > > > > After a little tweaking (the way AMBER names histidine residues), I > > successful converted the molecule (approx 3100 amino acids) into a CG > > representation. I successfully energy min the protein in vacuum to a > > threshold of 500, and in solvent to a threshold of 750 using steepest > > descent. Looking for a system at an energy min of a threshold around 300 > I > > begin to see LINCS warnings. Observing the initial structure, there is > > nothing obviously wrong with the bond network (both protein and polarized > > CG water). > > > > I take the system that energy mins at 750 (protein-water mix, with no > > fault reported), and went straight to NPT, 20fs step. Blew up. After a > bit > > of chatting with the MARTINI community, I¹ve started with an NVT > ensemble, > > beginning at 5s then through 10fs, 15fs, and 20fs. I only run for 1000 > > steps before switching. Keeping any of the simulations running for longer > > throws lincs warnings followed by a segmentation fault from the warning: > > > > "3 particles communicated to PME rank 7 are more than 2/3 times the > > cut-off out of the domain decomposition cell of their charge group in > > dimension x." > > > > Observing the trajectories of any of the extended simulations shows the > > protein snapping like a rope, and always at the same place. I have > watched > > every trajectory at this point, using numerous energy min start points, > to > > try and understand why it is blowing up. I can¹t see any obvious reason. > I > > was told to consider how the temperature is changing. Below is an example > > of the temperature and pressure from an NPT of 20fs step continued from > > the very short 20fs step NVT simulation (hoping that perhaps CG without > > pressure just doesn¹t behave happily; I was wrong). > > > > > > TEMP: > > Š > > 6.630000 311.000336 > > 6.645000 311.371643 > > 6.660000 311.724213 > > 6.675000 313.878693 > > 6.690000 681558.937500 > > > > > > PRESSURE: > > Š > > 6.630000 3.559879 > > 6.645000 3.901433 > > 6.660000 3.589078 > > 6.675000 4.158611 > > 6.690000 81762.437500 > > > > The final LINCS warning from this same run: > > > > Step 300, time 4.5 (ps) LINCS WARNING > > relative constraint deviation after LINCS: > > rms 0.000035, max 0.003386 (between atoms 2125 and 2126) > > bonds that rotated more than 45 degrees: > > atom 1 atom 2 angle previous, current, constraint length > > 2125 2126 68.3 0.2781 0.2691 0.2700 > > 2125 2127 45.9 0.2789 0.2701 0.2700 > > > > > > At this stage the structure ruptures as described above. > > > > > > My NVT settings (with NPT included to save space) are: > > > > ----------------- > > title = Martini > > > > integrator = md > > dt = 0.015 > > nsteps = 1000 > > nstcomm = 100 > > ;comm-grps = > > > > nstxout = 1000 > > nstvout = 1000 > > nstfout = 0 > > nstlog = 1 > > nstenergy = 1 > > nstxout-compressed = 0 > > compressed-x-precision = 0 > > ;compressed-x-grps = > > energygrps = collagen solvent > > > > cutoff-scheme = Verlet > > nstlist = 20 > > ns_type = grid > > > > pbc = xyz > > verlet-buffer-tolerance = 0.005 > > > > coulombtype = PME ;reaction-field > > rcoulomb = 1.1 > > fourierspacing = 0.16 ;0.2 ;0.12 > > > > epsilon_r = 2.5 ;15 ; 2.5 (with polarizable water) > > epsilon_rf = 0 > > vdw_type = cutoff > > vdw-modifier = Potential-shift-verlet > > rvdw = 1.1 > > > > tcoupl = v-rescale ;berendsen ;v-rescale > > tc-grps = collagen solvent > > tau_t = 0.5 0.5 ;1.0 1.0 > > ref_t = 310 310 > > > > Pcoupl = berendsen ;parrinello-rahman > > Pcoupltype = isotropic > > tau_p = 12.0 ; parrinello-rahman is more stable with > > larger tau-p, DdJ, 20130422 > > compressibility = 10e-4 > > ref_p = 1.0 > > refcoord_scaling = com > > > > gen_vel = no > > gen_temp = 310 > > gen_seed = 473529 > > > > continuation = yes > > constraints = none > > constraint_algorithm = lincs > > lincs-warnangle = 45 > > lincs-order=8 > > lincs-iter=4 > > > > > > ‹‹‹‹‹‹‹‹‹‹ > > > > Every setting bar the lincs iter, order, warnangle were supplied with the > > latest version of MARTINI. During many NVT runs I have adjusted the tau-t > > to try and keep the thermostat from oscillating its way into infinity. > > > > I¹m curious, will an out of control thermostat break a structure, or will > > a structure breaking (for what ever reason this structure is breaking) > > cause the thermostat to go out of control? > > > > My only thought thought is the initial .itp file that Martinize created. > I > > informed the script that this was collagen, therefor it sets ³F² and the > > corresponding bonded parameters. A human collagen is not perfect in its > > helical structure. Could there be underlying forces contributed from > badly > > bonded backbone-backbone arrangements? > > > > [ atoms ] > > 1 N0 1 GLY BB 1 0.0000 ; F > > 2 N0 2 LEU BB 2 0.0000 ; F > > 3 C1 2 LEU SC1 3 0.0000 ; F > > 4 N0 3 SER BB 4 0.0000 ; F > > > > > > Many thanks > > Anthony > > > > Thanks > > Anthony > > > > -- > > Gromacs Users mailing list > > > > * Please search the archive at > > http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before > > posting! > > > > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists > > > > * For (un)subscribe requests visit > > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or > > send a mail to gmx-users-requ...@gromacs.org. > > > -- > Gromacs Users mailing list > > * Please search the archive at http://www.gromacs.org/ > Support/Mailing_Lists/GMX-Users_List before posting! > > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists > > * For (un)subscribe requests visit > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or > send a mail to gmx-users-requ...@gromacs.org. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.