Re: [gmx-users] Lincs all-bonds causing instability in otherwise stable system
Hi, After further thought, the base problem is that there must be at least one rigid triangle, and LINCS can't do this because it couples constraints too much. The water models use [settles] to solve the particular case of a single triangle, and TIP4P and TIP5P water models are implemented with virtual sites on such a triangle. So with a bit of pencil and paper, I'm sure you could e.g. have a Me-Me-O rigid triangle and use a 3out virtual site (see manual) for the sulfur. But you'd only do this if the need to use SHAKE was hurting your simulation performance. Mark On Thu, May 21, 2015 at 10:16 AM Cara Kreck cara.kr...@student.curtin.edu.au wrote: Thanks Mark, A grompp warning about problematic constraints sounds like a good idea. DMSO is actually tetrahedral, not planar. Is that possible (and hopefully straightforward) to do with virtual sites? If it is, is there an example I could look at? Thanks, Cara From: mark.j.abra...@gmail.com Date: Wed, 20 May 2015 11:47:43 + To: gmx-us...@gromacs.org; gromacs.org_gmx-users@maillist.sys.kth.se Subject: Re: [gmx-users] Lincs all-bonds causing instability in otherwise stable system Hi, Sorry that was so painful for you. The LINCS documentation in the manual suggests that that combination of bonds becoming constraints won't work well. Perhaps we should add a 3-edge cycle detector to grompp? If you actually want a rigid planar DMSO, then three real atoms and a virtual site would be the best implementation. Mark On Wed, May 20, 2015 at 6:41 AM Cara Kreck cara.kr...@student.curtin.edu.au wrote: Hi everyone, I never got a reply to my message but I figured out the problem by myself. Just in case anyone else runs into a similar problem I thought I should explain what was wrong and share the solution. I was using a DMSO topology from the ATB that uses extra bonds to fix the geometry instead of angle or dihedral terms: [ moleculetype ] ; Name nrexcl DMSO 3 [ atoms ] ; nr type resnr resid atom cgnr chargemasstotal_charge 1 SDmso1DMSO S10.128 32.0600 2 ODmso1DMSO O1 -0.448 15.9994 3 CDmso1DMSO C10.160 15.0350 4 CDmso1DMSO C10.160 15.0350 ; 0.000 ; total charge of the molecule: 0.000 [ bonds ] ; ai aj funct c0 c1 122 0.1530 8.0400e+06 132 0.1938 4.9500e+06 142 0.1938 4.9500e+06 232 0.2794 2.3900e+06 242 0.2794 2.3900e+06 342 0.2912 2.1900e+06 This topology is fine to use with SHAKE but LINCS doesn't seem to be able to handle it. When I removed the extra bonds my simulations were able to run with all bonds constrained by LINCS. Then I found appropriate angle and dihedral parameters in the GROMOS ffbonded.itp file to control the geometry again. My topology file now looks like this: [ moleculetype ] ; Name nrexcl DMSO 3 [ atoms ] ; nr type resnr resid atom cgnr chargemasstotal_charge 1 SDmso1DMSO S10.128 32.0600 2 ODmso1DMSO O1 -0.448 15.9994 3 CDmso1DMSO C10.160 15.0350 4 CDmso1DMSO C10.160 15.0350 ; 0.000 ; total charge of the molecule: 0.000 [ bonds ] ; ai aj funct c0 c1 122 0.1530 8.0400e+06 132 0.1938 4.9500e+06 142 0.1938 4.9500e+06 [ angles ] ; ai aj ak funct angle fc 314297.40 469.00 3122 106.75 503.00 4122 106.75 503.00 [ dihedrals ] ; GROMOS improper dihedrals ; ai aj ak al funct angle fc 13422 35.26439 334.84617 ; tetrahedral centre I ran an energy minimisation of a single molecule with the new topology and its geometry overlapped the old one perfectly. So the problem was difficult to diagnose but easy to fix. Especially since I was able to energy minimise the individual molecules with all bonds constrained but the constraints went haywire when everything was combined in the full system. Hopefully I can at least save someone else from wasting 3 weeks trying to get a similar topology to work with LINCS. Cara From: cara.kr...@student.curtin.edu.au To: gromacs.org_gmx-users@maillist.sys.kth.se Date: Fri, 15 May 2015 14:31:16 +0800 Subject: [gmx-users] Lincs all-bonds causing instability in otherwise stable system I forgot to include an example of my mdp files. I've tried varying the timestep, running with and without pressure coupling, and obviously changing the constraints as outlined in the previous message
Re: [gmx-users] Lincs all-bonds causing instability in otherwise stable system
Thanks Mark, A grompp warning about problematic constraints sounds like a good idea. DMSO is actually tetrahedral, not planar. Is that possible (and hopefully straightforward) to do with virtual sites? If it is, is there an example I could look at? Thanks, Cara From: mark.j.abra...@gmail.com Date: Wed, 20 May 2015 11:47:43 + To: gmx-us...@gromacs.org; gromacs.org_gmx-users@maillist.sys.kth.se Subject: Re: [gmx-users] Lincs all-bonds causing instability in otherwise stable system Hi, Sorry that was so painful for you. The LINCS documentation in the manual suggests that that combination of bonds becoming constraints won't work well. Perhaps we should add a 3-edge cycle detector to grompp? If you actually want a rigid planar DMSO, then three real atoms and a virtual site would be the best implementation. Mark On Wed, May 20, 2015 at 6:41 AM Cara Kreck cara.kr...@student.curtin.edu.au wrote: Hi everyone, I never got a reply to my message but I figured out the problem by myself. Just in case anyone else runs into a similar problem I thought I should explain what was wrong and share the solution. I was using a DMSO topology from the ATB that uses extra bonds to fix the geometry instead of angle or dihedral terms: [ moleculetype ] ; Name nrexcl DMSO 3 [ atoms ] ; nr type resnr resid atom cgnr chargemasstotal_charge 1 SDmso1DMSO S10.128 32.0600 2 ODmso1DMSO O1 -0.448 15.9994 3 CDmso1DMSO C10.160 15.0350 4 CDmso1DMSO C10.160 15.0350 ; 0.000 ; total charge of the molecule: 0.000 [ bonds ] ; ai aj funct c0 c1 122 0.1530 8.0400e+06 132 0.1938 4.9500e+06 142 0.1938 4.9500e+06 232 0.2794 2.3900e+06 242 0.2794 2.3900e+06 342 0.2912 2.1900e+06 This topology is fine to use with SHAKE but LINCS doesn't seem to be able to handle it. When I removed the extra bonds my simulations were able to run with all bonds constrained by LINCS. Then I found appropriate angle and dihedral parameters in the GROMOS ffbonded.itp file to control the geometry again. My topology file now looks like this: [ moleculetype ] ; Name nrexcl DMSO 3 [ atoms ] ; nr type resnr resid atom cgnr chargemasstotal_charge 1 SDmso1DMSO S10.128 32.0600 2 ODmso1DMSO O1 -0.448 15.9994 3 CDmso1DMSO C10.160 15.0350 4 CDmso1DMSO C10.160 15.0350 ; 0.000 ; total charge of the molecule: 0.000 [ bonds ] ; ai aj funct c0 c1 122 0.1530 8.0400e+06 132 0.1938 4.9500e+06 142 0.1938 4.9500e+06 [ angles ] ; ai aj ak funct angle fc 314297.40 469.00 3122 106.75 503.00 4122 106.75 503.00 [ dihedrals ] ; GROMOS improper dihedrals ; ai aj ak al funct angle fc 13422 35.26439 334.84617 ; tetrahedral centre I ran an energy minimisation of a single molecule with the new topology and its geometry overlapped the old one perfectly. So the problem was difficult to diagnose but easy to fix. Especially since I was able to energy minimise the individual molecules with all bonds constrained but the constraints went haywire when everything was combined in the full system. Hopefully I can at least save someone else from wasting 3 weeks trying to get a similar topology to work with LINCS. Cara From: cara.kr...@student.curtin.edu.au To: gromacs.org_gmx-users@maillist.sys.kth.se Date: Fri, 15 May 2015 14:31:16 +0800 Subject: [gmx-users] Lincs all-bonds causing instability in otherwise stable system I forgot to include an example of my mdp files. I've tried varying the timestep, running with and without pressure coupling, and obviously changing the constraints as outlined in the previous message: integrator = md dt = 0.001 ; 1fs nsteps = 10 ; 100ps comm_grps= DOPC !DOPC nstxout = 1000 nstvout = 0 nstlog = 1000 nstenergy= 1000 energygrps = DOPC !DOPC nstcalcenergy= 5 nstlist = 5 ns_type = grid pbc = xyz rlist= 0.8 coulombtype = Reaction-Field rcoulomb = 1.4 epsilon_rf = 62 vdwtype = Cut-off rvdw = 1.4 tcoupl = berendsen tc-grps = DOPC !DOPC
Re: [gmx-users] Lincs all-bonds causing instability in otherwise stable system
Hi, Sorry that was so painful for you. The LINCS documentation in the manual suggests that that combination of bonds becoming constraints won't work well. Perhaps we should add a 3-edge cycle detector to grompp? If you actually want a rigid planar DMSO, then three real atoms and a virtual site would be the best implementation. Mark On Wed, May 20, 2015 at 6:41 AM Cara Kreck cara.kr...@student.curtin.edu.au wrote: Hi everyone, I never got a reply to my message but I figured out the problem by myself. Just in case anyone else runs into a similar problem I thought I should explain what was wrong and share the solution. I was using a DMSO topology from the ATB that uses extra bonds to fix the geometry instead of angle or dihedral terms: [ moleculetype ] ; Name nrexcl DMSO 3 [ atoms ] ; nr type resnr resid atom cgnr chargemasstotal_charge 1 SDmso1DMSO S10.128 32.0600 2 ODmso1DMSO O1 -0.448 15.9994 3 CDmso1DMSO C10.160 15.0350 4 CDmso1DMSO C10.160 15.0350 ; 0.000 ; total charge of the molecule: 0.000 [ bonds ] ; ai aj funct c0 c1 122 0.1530 8.0400e+06 132 0.1938 4.9500e+06 142 0.1938 4.9500e+06 232 0.2794 2.3900e+06 242 0.2794 2.3900e+06 342 0.2912 2.1900e+06 This topology is fine to use with SHAKE but LINCS doesn't seem to be able to handle it. When I removed the extra bonds my simulations were able to run with all bonds constrained by LINCS. Then I found appropriate angle and dihedral parameters in the GROMOS ffbonded.itp file to control the geometry again. My topology file now looks like this: [ moleculetype ] ; Name nrexcl DMSO 3 [ atoms ] ; nr type resnr resid atom cgnr chargemasstotal_charge 1 SDmso1DMSO S10.128 32.0600 2 ODmso1DMSO O1 -0.448 15.9994 3 CDmso1DMSO C10.160 15.0350 4 CDmso1DMSO C10.160 15.0350 ; 0.000 ; total charge of the molecule: 0.000 [ bonds ] ; ai aj funct c0 c1 122 0.1530 8.0400e+06 132 0.1938 4.9500e+06 142 0.1938 4.9500e+06 [ angles ] ; ai aj ak funct angle fc 314297.40 469.00 3122 106.75 503.00 4122 106.75 503.00 [ dihedrals ] ; GROMOS improper dihedrals ; ai aj ak al funct angle fc 13422 35.26439 334.84617 ; tetrahedral centre I ran an energy minimisation of a single molecule with the new topology and its geometry overlapped the old one perfectly. So the problem was difficult to diagnose but easy to fix. Especially since I was able to energy minimise the individual molecules with all bonds constrained but the constraints went haywire when everything was combined in the full system. Hopefully I can at least save someone else from wasting 3 weeks trying to get a similar topology to work with LINCS. Cara From: cara.kr...@student.curtin.edu.au To: gromacs.org_gmx-users@maillist.sys.kth.se Date: Fri, 15 May 2015 14:31:16 +0800 Subject: [gmx-users] Lincs all-bonds causing instability in otherwise stable system I forgot to include an example of my mdp files. I've tried varying the timestep, running with and without pressure coupling, and obviously changing the constraints as outlined in the previous message: integrator = md dt = 0.001 ; 1fs nsteps = 10 ; 100ps comm_grps= DOPC !DOPC nstxout = 1000 nstvout = 0 nstlog = 1000 nstenergy= 1000 energygrps = DOPC !DOPC nstcalcenergy= 5 nstlist = 5 ns_type = grid pbc = xyz rlist= 0.8 coulombtype = Reaction-Field rcoulomb = 1.4 epsilon_rf = 62 vdwtype = Cut-off rvdw = 1.4 tcoupl = berendsen tc-grps = DOPC !DOPC tau_t= 0.1 0.1 ref_t= 303 303 ;Pcoupl = berendsen ;pcoupltype = semiisotropic ;tau_p= 1.0 1.0 ;compressibility = 4.6e-5 4.6e-5 ;ref_p= 1.0 1.0 gen_vel = no constraints = all-bonds constraint_algorithm = shake continuation = yes From: cara.kr...@student.curtin.edu.au To: gromacs.org_gmx-users@maillist.sys.kth.se Date: Fri, 15 May 2015 14:17:28 +0800 Subject: [gmx-users] FW: Lincs all-bonds causing instability in
Re: [gmx-users] Lincs all-bonds causing instability in otherwise stable system
Hi everyone, I never got a reply to my message but I figured out the problem by myself. Just in case anyone else runs into a similar problem I thought I should explain what was wrong and share the solution. I was using a DMSO topology from the ATB that uses extra bonds to fix the geometry instead of angle or dihedral terms: [ moleculetype ] ; Name nrexcl DMSO 3 [ atoms ] ; nr type resnr resid atom cgnr chargemasstotal_charge 1 SDmso1DMSO S10.128 32.0600 2 ODmso1DMSO O1 -0.448 15.9994 3 CDmso1DMSO C10.160 15.0350 4 CDmso1DMSO C10.160 15.0350 ; 0.000 ; total charge of the molecule: 0.000 [ bonds ] ; ai aj funct c0 c1 122 0.1530 8.0400e+06 132 0.1938 4.9500e+06 142 0.1938 4.9500e+06 232 0.2794 2.3900e+06 242 0.2794 2.3900e+06 342 0.2912 2.1900e+06 This topology is fine to use with SHAKE but LINCS doesn't seem to be able to handle it. When I removed the extra bonds my simulations were able to run with all bonds constrained by LINCS. Then I found appropriate angle and dihedral parameters in the GROMOS ffbonded.itp file to control the geometry again. My topology file now looks like this: [ moleculetype ] ; Name nrexcl DMSO 3 [ atoms ] ; nr type resnr resid atom cgnr chargemasstotal_charge 1 SDmso1DMSO S10.128 32.0600 2 ODmso1DMSO O1 -0.448 15.9994 3 CDmso1DMSO C10.160 15.0350 4 CDmso1DMSO C10.160 15.0350 ; 0.000 ; total charge of the molecule: 0.000 [ bonds ] ; ai aj funct c0 c1 122 0.1530 8.0400e+06 132 0.1938 4.9500e+06 142 0.1938 4.9500e+06 [ angles ] ; ai aj ak funct angle fc 314297.40 469.00 3122 106.75 503.00 4122 106.75 503.00 [ dihedrals ] ; GROMOS improper dihedrals ; ai aj ak al funct angle fc 13422 35.26439 334.84617 ; tetrahedral centre I ran an energy minimisation of a single molecule with the new topology and its geometry overlapped the old one perfectly. So the problem was difficult to diagnose but easy to fix. Especially since I was able to energy minimise the individual molecules with all bonds constrained but the constraints went haywire when everything was combined in the full system. Hopefully I can at least save someone else from wasting 3 weeks trying to get a similar topology to work with LINCS. Cara From: cara.kr...@student.curtin.edu.au To: gromacs.org_gmx-users@maillist.sys.kth.se Date: Fri, 15 May 2015 14:31:16 +0800 Subject: [gmx-users] Lincs all-bonds causing instability in otherwise stable system I forgot to include an example of my mdp files. I've tried varying the timestep, running with and without pressure coupling, and obviously changing the constraints as outlined in the previous message: integrator = md dt = 0.001 ; 1fs nsteps = 10 ; 100ps comm_grps= DOPC !DOPC nstxout = 1000 nstvout = 0 nstlog = 1000 nstenergy= 1000 energygrps = DOPC !DOPC nstcalcenergy= 5 nstlist = 5 ns_type = grid pbc = xyz rlist= 0.8 coulombtype = Reaction-Field rcoulomb = 1.4 epsilon_rf = 62 vdwtype = Cut-off rvdw = 1.4 tcoupl = berendsen tc-grps = DOPC !DOPC tau_t= 0.1 0.1 ref_t= 303 303 ;Pcoupl = berendsen ;pcoupltype = semiisotropic ;tau_p= 1.0 1.0 ;compressibility = 4.6e-5 4.6e-5 ;ref_p= 1.0 1.0 gen_vel = no constraints = all-bonds constraint_algorithm = shake continuation = yes From: cara.kr...@student.curtin.edu.au To: gromacs.org_gmx-users@maillist.sys.kth.se Date: Fri, 15 May 2015 14:17:28 +0800 Subject: [gmx-users] FW: Lincs all-bonds causing instability in otherwise stable system Hi, I am having trouble constraining all-bonds with lincs in a system that is stable using shake all-bonds or lincs h-bonds in gromacs 4.6.7. The previous step using lincs h-bonds gave these energies: Step Time Lambda 20 200.00.0 Energies (kJ/mol) G96Bond G96AngleProper Dih. Improper Dih. LJ-14 1.76317e+042.43132e+041.39308e+04
