[gmx-users] Cyclic Peptide in Cyclohexane - LINCS warnings
Billy Williams-Noonan
billy.williams-noonan at monash.edu
Mon Jun 12 10:59:35 CEST 2017
Hi All,
I've increased tau_p to 6 ps with a 2 fs time-step, and the systems do not
crash / generate LINCS warnings. My background is not in mathematics. Is
this solution a reasonable thing to do (at least more so than using the
Berendsen barostat)?
Cheers,
Billy
On 12 June 2017 at 13:41, Billy Williams-Noonan <
billy.williams-noonan at monash.edu> wrote:
> Hi GROMACS Users,
>
> Am running a simulation of a cyclic peptide with some unusual residues in
> a box of cyclohexane.
>
> The parameters for this peptide are sources from the ATB, as are the
> parameters for the cyclohexane solvent. Am using the Gromos 54a7
> parameters with GROMACS-2016.
>
> My simulation keeps crashing after repeated LINCS warnings and after
> re-running the simulations with just a box of cyclohexane, as well as
> re-running with a simulation of the cyclic peptide in vacuum, I've found
> that it is the cyclohexane, specifically, that is the problem. By looking
> at the step files, you can see that the atoms of each solvent molecule come
> into close-contact with each other, creating a high amount of force to be
> experienced by the relevant atoms, and thereby causing the relevant
> molecules to 'blow up'.
>
> I was wondering why this would be the case.
>
> Here are my *.mdp *details.
>
> title = test_cylohexane
> ; Run parameters
> integrator = md ; leap-frog integrator
> nsteps = 50000000 ; = 2 * 50000 ps = 100 ns
> dt = 0.002 ; 1 fs
> ; Output control
> nstxout = 5000 ; save coordinates every 1.0 ps
> nstvout = 5000 ; save velocities every 1.0 ps
> nstenergy = 5000 ; save energies every 1.0 ps
> nstlog = 5000 ; update log file every 1.0 ps
> ; Bond parameters
> continuation = yes ; Restarting after NVT
> constraint_algorithm = lincs ; holonomic constraints
> constraints = all-bonds ; all bonds (even heavy atom-H
> bonds) constrained
> lincs_iter = 1 ; accuracy of LINCS
> lincs_order = 4 ; also related to accuracy
> ; Neighborsearching
> cutoff-scheme = Verlet
> ns_type = grid ; search neighboring grid cells
> nstlist = 40 ; 80 fs, largely irrelevant with
> Verlet scheme
> rcoulomb = 1.0 ; short-range electrostatic cutoff
> (in nm)
> rvdw = 1.0 ; short-range van der Waals
> cutoff (in nm)
> ; Electrostatics
> coulombtype = PME ; Particle Mesh Ewald for
> long-range electrostatics
> pme_order = 4 ; cubic interpolation
> fourierspacing = 0.16 ; grid spacing for FFT
> ; Temperature coupling is on
> tcoupl = V-rescale ; modified Berendsen
> thermostat
> tc-grps = half_A half_B ; two coupling groups - more
> accurate
> tau_t = 0.1 0.1 ; time constant, in ps
> ref_t = 300 300 ; reference temperature, one for
> each group, in K
> ; Pressure coupling is on
> pcoupl = Parrinello-Rahman ; Pressure coupling on
> in NPT
> pcoupltype = isotropic ; uniform scaling of
> box vectors
> tau_p = 2.0 ; time constant, in ps
> ref_p = 1.0 ; reference pressure,
> in bar
> compressibility = 4.5e-5 ; isothermal compressibility
> of water, bar^-1
> ; Periodic boundary conditions
> pbc = xyz ; 3-D PBC
> ; Dispersion correction
> DispCorr = EnerPres ; account for cut-off vdW scheme
> ; Velocity generation
> gen_vel = no ; Velocity generation is off
>
> Here are is the topology for cyclohexane.
>
> [ moleculetype ]
> ; Name nrexcl
> _VMQ 3
> [ atoms ]
> ; nr type resnr resid atom cgnr charge mass total_charge
> 1 CH2r 1 _VMQ C1 1 0.000 14.0270
> 2 CH2r 1 _VMQ C6 1 0.000 14.0270 ; 0.000
> 3 CH2r 1 _VMQ C2 2 0.000 14.0270
> 4 CH2r 1 _VMQ C3 2 0.000 14.0270 ; 0.000
> 5 CH2r 1 _VMQ C4 3 0.000 14.0270
> 6 CH2r 1 _VMQ C5 3 0.000 14.0270 ; 0.000
> ; total charge of the molecule: 0.000
> [ bonds ]
> ; ai aj funct c0 c1
> 1 2 2 0.1530 7.1500e+06
> 1 3 2 0.1530 7.1500e+06
> 2 6 2 0.1530 7.1500e+06
> 3 4 2 0.1530 7.1500e+06
> 4 5 2 0.1530 7.1500e+06
> 5 6 2 0.1530 7.1500e+06
> [ pairs ]
> ; ai aj funct ; all 1-4 pairs but the ones excluded in GROMOS itp
> 1 5 1
> 2 4 1
> 3 6 1
> [ angles ]
> ; ai aj ak funct angle fc
> 2 1 3 2 111.00 530.00
> 1 2 6 2 111.00 530.00
> 1 3 4 2 111.00 530.00
> 3 4 5 2 111.00 530.00
> 4 5 6 2 111.00 530.00
> 2 6 5 2 111.00 530.00
> [ dihedrals ]
> ; GROMOS improper dihedrals
> ; ai aj ak al funct angle fc
> [ dihedrals ]
> ; ai aj ak al funct ph0 cp mult
> 3 1 2 6 1 0.00 5.92 3
> 2 1 3 4 1 0.00 5.92 3
> 1 2 6 5 1 0.00 5.92 3
> 1 3 4 5 1 0.00 5.92 3
> 3 4 5 6 1 0.00 5.92 3
> 4 5 6 2 1 0.00 5.92 3
> [ exclusions ]
> ; ai aj funct ; GROMOS 1-4 exclusions
>
>
> The periodic box dimensions were set up so that the density of the system
> was 0.779 g/mL, which is the density of cyclohexane from experiment.
>
> I should probably add that there was one simulation where I could not see
> two cyclohexane molecules obviously 'touching' (this was run in triplicate
> for 100 ns) - but two molecules of cyclohexane 'blew up' regardless.
> Further, reducing the time-step to 1 fs, or using a Berendsen barostat
> eliminates the LINCS warnings an prevents crashing, but I'm unsure if I'm
> actually fixing the problem by doing that. Further, my understanding is
> that Berendsen gives the wrong ensemble pressure, and I'd rather avoid
> using it, just as I'd like my simulations to run more quickly and have a 2
> fs time-step.
>
> Is there anything wrong with what I'm doing or anything I should add to
> the *.mdp* file to make it work?
>
> Cheers,
>
> Billy
>
>
>
>
> --
> Billy Noonan* | *PhD Student *|* Bsci ( *Adv* ), IA Hon
>
> *LinkedIn Profile
> <http://www.linkedin.com/profile/preview?locale=en_US&trk=prof-0-sb-preview-primary-button>
> **|* +61420 382 557 <+61%20420%20382%20557>
>
> Monash Institute for Pharmaceutical Sciences ( *MIPS* )
> Royal Parade, Parkville, 3052
>
>
--
Billy Noonan* | *PhD Student *|* Bsci ( *Adv* ), IA Hon
*LinkedIn Profile
<http://www.linkedin.com/profile/preview?locale=en_US&trk=prof-0-sb-preview-primary-button>
**|* +61420 382 557
Monash Institute for Pharmaceutical Sciences ( *MIPS* )
Royal Parade, Parkville, 3052
More information about the gromacs.org_gmx-users
mailing list