Could someone please check if the following minimization input file with harmonic restraints and the corresponding DISANG file is constructed correctly? I want to treat the Zn-coordinated water and glutamate with the restrained non-bonded model and the Zn-coordinated histidines with the bonded model
min.in:
Constant Volume Minimization 1 of solvent
# Control section
&cntrl
imin=1,
irest = 0,
ntmin = 1, maxcyc = 5000, ncyc = 2500, dx0 = 0.01, drms = 0.0001,
cut = 10.0,
ntb = 1,
ntp = 0,
ntwx = 50,
ntpr = 50,
ntwr = 5000,
ioutfm = 1,
iwrap = 1,
ntr = 1, restraintmask='!:WAT,Na+', restraint_wt=50.0,
nmropt = 1,
&end
&wt
type='END',
&end
DISANG=NDOM_rnb.RST
/
NDOM_rnb.RST:
# Harmonic restraints for the restrained nonbonded model of ACE NDOM
&rst iat=6143,12493, r1=0., r2=1.959, r3=1.959, r4=100., rk2=105.73, rk3=105.73,/
&rst iat=12528,12493, r1=0., r2=2.056, r3=2.056, r4=100., rk2=66.54, rk3=66.54,/
&end
Any advice would be much appreciated
Regards
Lizelle Lubbe
PhD (Medical biochemistry) candidate
Department of Integrative Biomedical Sciences
University of Cape Town
________________________________________
From: Lizelle Lubbe
Sent: 26 September 2017 12:34 PM
To: Pengfei Li
Subject: Re: [AMBER] SHAKE failure
Hi Pengfei,
Thanks for the reply last week. What you suggested makes sense and I want to try it out. I've read through the tutorials and the manual's section on NMR restraints but am unsure about a few things.
I think it might work if I apply the restrained non-bonded model to the glutamate and water molecule and then treat the two histidines still with the bonded model. Do you agree? Would I then just remove the bond information of these two Zn-bonds from the leap.in file and edit the tutorial's RST file to contain my atom numbers and force-constants?
As far as I understand I should copy the force-constants and bond lengths from the frcmod file used in the bonded model into the rk2=rk3 and r2=r3 columns, respectively. Is this correct? Should I set r1=0 and r4=100.0 as in the tutorial or modify them?
You mention in the tutorial that these harmonic restraints should be applied during minimization and dynamics. Could you please provide me with an example of how to source this RST file or include this information in the minimization.in file? I can't find anything in the manual...
Kind regards
Lizelle
________________________________________
From: Pengfei Li <ambermailpengfei.gmail.com>
Sent: 23 September 2017 11:34:31 PM
To: AMBER Mailing List
Subject: Re: [AMBER] SHAKE failure
Hi Lizelle,
I guess you met the problem mentioned in the end of the webpage: http://ambermd.org/tutorials/advanced/tutorial20/mcpbpy.htm <http://ambermd.org/tutorials/advanced/tutorial20/mcpbpy.htm>. To keep it short, this is because that the hydrogen atoms in water has zero VDW interaction with surroundings. However, since they have positive charges, when they have a 1-4 interaction any atoms with negative charges, the structural collapse would happen. This is because there is only 1-4 electrostatic attraction but no 1-4 VDW repulsion under this situation, this is an error can be meet in the bonded model of metal ion. To solve the problem, you can consider to use the restrained nonbonded model: http://ambermd.org/tutorials/advanced/tutorial20/rnb.htm <http://ambermd.org/tutorials/advanced/tutorial20/rnb.htm>.
Kind regards,
Pengfei
> On Sep 22, 2017, at 1:40 PM, Lizelle Lubbe <LBBLIZ002.myuct.ac.za> wrote:
>
> Sorry just a slight correction... During minimization 2 restraintmask :WAR & !.H=
>
> Lizelle
> ________________________________________
> From: Lizelle Lubbe <LBBLIZ002.myuct.ac.za>
> Sent: 22 September 2017 08:37:05 PM
> To: amber.ambermd.org; Daniel Roe; ross.cgl.ucsf.edu
> Subject: Re: [AMBER] SHAKE failure
>
> Dear Daniel and Bill,
>
> Thanks so much for your replies. Daniel, I have repeated the equilibration with printing at each step and saw that just as in the minimizations when having the water unrestrained, the oxygen moves sideways and causes an overlap of the H with OE1 of the Zn-bound glutamate. The initial pdb input for leap and leap output look fine and cpptraj checks out with no overlap messages. The overlap only occurs when allowing the water to move freely. I suspect that there is something wrong in my water parameters somewhere but don't know where to start searching. I've attached some images of the equilibration printing every step and some energy output below.
> Bill, the input files of all steps are given below. Also the lib file of the water molecule. I obtained this from another group and perhaps don't understand the content very well. I replaced the coordinates under !entry.'resname'.unit.positions table with my Zn-bound water coordinates from the pdb. Should I rather have retained the original values matching the ones for WAT in solvents.lib? I also replaced the Zn-bound histidine and glutamate coordinates with my own and am wondering if this was a mistake and results in this strange behaviour?
>
> original water lib file from publisher's of Zn parameters:
>
> !!index array str
> "WAR"
> !entry.WAR.unit.atoms table str name str type int typex int resx int flags int seq int elmnt dbl chg
> "O" "OK" 0 1 131072 3 8 -0.751392
> "H1" "HW" 0 1 131072 1 1 0.396784
> "H2" "HW" 0 1 131072 2 1 0.424564
> !entry.WAR.unit.atomspertinfo table str pname str ptype int ptypex int pelmnt dbl pchg
> "O" "OK" 0 -1 0.0
> "H1" "HW" 0 -1 0.0
> "H2" "HW" 0 -1 0.0
> !entry.WAR.unit.boundbox array dbl
> -1.000000
> 0.0
> 0.0
> 0.0
> 0.0
> !entry.WAR.unit.childsequence single int
> 2
> !entry.WAR.unit.connect array int
> 0
> 0
> !entry.WAR.unit.connectivity table int atom1x int atom2x int flags
> 3 1 1
> 2 1 1
> 2 3 1
> !entry.WAR.unit.hierarchy table str abovetype int abovex str belowtype int belowx
> "U" 0 "R" 1
> "R" 1 "A" 1
> "R" 1 "A" 3
> "R" 1 "A" 2
> !entry.WAR.unit.name single str
> "TP3"
> !entry.WAR.unit.positions table dbl x dbl y dbl z
> 0.0 0.0 0.0
> 0.957200 0.0 0.0
> -0.239988 0.926627 0.0
> !entry.WAR.unit.residueconnect table int c1x int c2x int c3x int c4x int c5x int c6x
> 0 0 0 0 0 0
> !entry.WAR.unit.residues table str name int seq int childseq int startatomx str restype int imagingx
> "WAR" 1 4 1 "w" 1
> !entry.WAR.unit.residuesPdbSequenceNumber array int
> 0
> !entry.WAR.unit.solventcap array dbl
> -1.000000
> 0.0
> 0.0
> 0.0
> 0.0
> !entry.WAR.unit.velocities table dbl x dbl y dbl z
> 0.0 0.0 0.0
> 0.0 0.0 0.0
> 0.0 0.0 0.0
>
> minimization 1 of solvent:
>
> # Control section
> &cntrl
> imin=1,
> irest = 0,
> ntmin = 1, maxcyc = 5000, ncyc = 2500, dx0 = 0.01, drms = 0.0001,
> cut = 10.0,
> ntb = 1,
> ntp = 0,
> ntwx = 50,
> ntpr = 50,
> ntwr = 5000,
> ioutfm = 1,
> iwrap = 1,
> ntr = 1, restraintmask='!:WAT,Na+', restraint_wt=50.0
>
> /
>
> minimization 2 of system (WAR = Zn-bound water):
>
> # Control section
> &cntrl
> imin=1,
> irest = 0,
> ntmin = 1, maxcyc = 20000, ncyc = 10000, dx0 = 0.01, drms = 0.0001,
> cut = 10.0,
> ntb = 1,
> ntp = 0,
> ntwx = 100,
> ntpr = 100,
> ntwr = 20000,
> ioutfm = 1,
> iwrap = 1,
> ntr = 1, restraintmask=':WAR', restraint_wt=50.0
>
> /
>
> heating solvent:
>
> # Control section
> &cntrl
> imin = 0,
> irest = 0,
> ntx = 1,
> nstlim = 300000, dt = 0.002,
> ntc = 2, ntf = 2,
> tempi = 10.0, temp0 = 300.0,
> ntwx = 5000, ntpr = 5000, ntwe = 5000,
> ntb = 1, ntp = 0,
> ntt = 3, gamma_ln = 2, ig=-1,
> cut = 10.0,
> ntxo = 2,
> ioutfm = 1,
> iwrap = 1,
> nmropt = 1,
> ntr = 1, restraintmask='!:WAT,Na+', restraint_wt=50.0
> /
> &wt type='TEMP0', istep1=0, istep2=100000, value1=10.0, value2=100.0/
> &wt type='TEMP0', istep1=100001, istep2=200000, value1=100.0, value2=200.0/
> &wt type='TEMP0', istep1=200001, istep2=300000, value1=200.0, value2=300.0/
> &wt type='END'/
>
> equilibration 1 of solvent:
>
> # Control section
> &cntrl
> imin = 0,
> irest = 1,
> ntx = 5,
> ntb = 2, ntp = 1, pres0 = 1.0, taup = 2.0,
> cut = 10.0,
> ntc = 2, ntf = 2,
> temp0 = 300.0,
> ntt = 3, gamma_ln = 2, ig = -1,
> nstlim = 100000, dt = 0.002,
> ntwx = 1000, ntpr = 1000, ntwe = 1000,
> ntwr = 100000,
> ioutfm=1,
> iwrap=1,
> ntr = 1, restraintmask='!:WAT,Na+', restraint_wt=10.0,
> /
>
> equilibration 2 of system:
>
> # Control section
> &cntrl
> imin = 0,
> irest = 1,
> ntx = 5,
> ntb = 2, ntp = 1, pres0 = 1.0, taup = 2.0,
> cut = 10.0,
> ntc = 2, ntf = 2,
> temp0 = 300.0,
> ntt = 3, gamma_ln = 2, ig = -1,
> nstlim = 100, dt = 0.002,
> ntwx = 1, ntpr = 1, ntwe = 1,
> ntwr = 100,
> ioutfm=1,
> iwrap=1,
> ntr = 0,
> /
>
> Selected output of eq2:
>
> NSTEP = 1 TIME(PS) = 800.002 TEMP(K) = 299.47 PRESS = 111.4
> Etot = -457615.5261 EKtot = 115240.5321 EPtot = -572856.0582
> BOND = 2615.1688 ANGLE = 6117.2451 DIHED = 7678.7143
> 1-4 NB = 3205.5098 1-4 EEL = 46188.6133 VDWAALS = 79065.9051
> EELEC = -717727.2146 EHBOND = 0.0000 RESTRAINT = 0.0000
> EKCMT = 53294.0126 VIRIAL = 48629.1977 VOLUME = 1939919.9138
> Density = 0.9935
> Ewald error estimate: 0.4355E-04
> ------------------------------------------------------------------------------
> NSTEP = 20 TIME(PS) = 800.040 TEMP(K) = 297.96 PRESS = 23.7
> Etot = -457190.3038 EKtot = 114659.3744 EPtot = -571849.6782
> BOND = 2699.5720 ANGLE = 6578.6660 DIHED = 7858.9118
> 1-4 NB = 3248.6770 1-4 EEL = 46041.9874 VDWAALS = 79365.4920
> EELEC = -717642.9844 EHBOND = 0.0000 RESTRAINT = 0.0000
> EKCMT = 53460.8164 VIRIAL = 52469.7550 VOLUME = 1939949.3169
> Density = 0.9935
> Ewald error estimate: 0.5955E-04
> ------------------------------------------------------------------------------
> NSTEP = 40 TIME(PS) = 800.080 TEMP(K) = 299.86 PRESS = -9.9
> Etot = -455438.8761 EKtot = 115389.4418 EPtot = -570828.3179
> BOND = 3420.4131 ANGLE = 6819.4033 DIHED = 7909.4072
> 1-4 NB = 3488.2059 1-4 EEL = 45783.8401 VDWAALS = 78987.9064
> EELEC = -717237.4939 EHBOND = 0.0000 RESTRAINT = 0.0000
> EKCMT = 53403.0882 VIRIAL = 53819.6086 VOLUME = 1939986.1076
> Density = 0.9935
> Ewald error estimate: 0.4069E-04
> ------------------------------------------------------------------------------
>
> vlimit exceeded for step 40; vmax = 22.4968
>
> Coordinate resetting (SHAKE) cannot be accomplished,
> deviation is too large
> NITER, NIT, LL, I and J are : 0 1 6057 12494 12496
>
> Thanks for the assistance
>
> Lizelle Lubbe
>
> PhD (Medical biochemistry) candidate
> Department of Integrative Biomedical Sciences
> University of Cape Town
> ________________________________________
> From: Bill Ross <ross.cgl.ucsf.edu>
> Sent: 22 September 2017 04:18 PM
> To: amber.ambermd.org
> Subject: Re: [AMBER] SHAKE failure
>
> Please copy/paste your min.in and md.in files for minimization and heating.
>
> Bill
>
>
> On 9/22/17 7:08 AM, Daniel Roe wrote:
>> On Fri, Sep 22, 2017 at 8:30 AM, Lizelle Lubbe <LBBLIZ002.myuct.ac.za> wrote:
>>> I've prepared my glycosylated metalloprotein using tleap and attempted minimization, heating and equilibration. Upon equilibration (having ntc = ntf = 2) I almost instantly get the following error:
>>>
>>> vlimit exceeded for step 61; vmax = 27.0336
>>> vlimit exceeded for step 62; vmax = 536.4434
>>>
>>> Coordinate resetting (SHAKE) cannot be accomplished,
>>> deviation is too large
>>> NITER, NIT, LL, I and J are : 0 5 6058 12494 12495
>>>
>>> Note: This is usually a symptom of some deeper
>>> problem with the energetics of the system.
>>>
>>> Atom 12494 is the O and atom 12495 the H1 of a Zn-coordinated water molecule (described using published forcefield parameters). Since the O and H1 are explicitly bonded I do not understand why SHAKE would want to constrain these bonds. Am I reading the error correctly? In parmed.py printBonds, printAngles, printDetails and printDihedrals for this water molecule all seem fine.
>> When you specify ntc = 2, *all* bonds to hydrogen are constrained. The
>> vlimit errors indicate your system is blowing up. It may be worthwhile
>> to re-run the MD and print out energies and trajectory frames every
>> step so you can see what is happening in more detail (especially
>> because it happens so quickly the file sizes should remain
>> reasonable). What may be happening is that the H of the water ends up
>> having a bad overlap with the Zn and shoots off into oblivion or
>> something. One thing you could try to go is check for bad overlaps in
>> your initial structure using the 'check' command from cpptraj:
>> https://github.com/Amber-MD/cpptraj. You may need more minimization
>> before running your MD phase.
>>
>> -Dan
>>
>>> Could anyone please help me to solve this issue?
>>>
>>> Lizelle Lubbe
>>>
>>> PhD (Medical biochemistry) candidate
>>> Department of Integrative Biomedical Sciences
>>> University of Cape Town
>>> Disclaimer - University of Cape Town This e-mail is subject to UCT policies and e-mail disclaimer published on our website at http://www.uct.ac.za/about/policies/emaildisclaimer/ or obtainable from +27 21 650 9111. If this e-mail is not related to the business of UCT, it is sent by the sender in an individual capacity. Please report security incidents or abuse via csirt.uct.ac.za
>>>
>>> _______________________________________________
>>> AMBER mailing list
>>> AMBER.ambermd.org
>>> http://lists.ambermd.org/mailman/listinfo/amber
>>
>>
>
>
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
> Disclaimer - University of Cape Town This e-mail is subject to UCT policies and e-mail disclaimer published on our website at http://www.uct.ac.za/about/policies/emaildisclaimer/ or obtainable from +27 21 650 9111. If this e-mail is not related to the business of UCT, it is sent by the sender in an individual capacity. Please report security incidents or abuse via csirt.uct.ac.za
> Disclaimer - University of Cape Town This e-mail is subject to UCT policies and e-mail disclaimer published on our website at http://www.uct.ac.za/about/policies/emaildisclaimer/ or obtainable from +27 21 650 9111. If this e-mail is not related to the business of UCT, it is sent by the sender in an individual capacity. Please report security incidents or abuse via csirt.uct.ac.za
>
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Disclaimer - University of Cape Town This e-mail is subject to UCT policies and e-mail disclaimer published on our website at
http://www.uct.ac.za/about/policies/emaildisclaimer/ or obtainable from +27 21 650 9111. If this e-mail is not related to the business of UCT, it is sent by the sender in an individual capacity. Please report security incidents or abuse via csirt.uct.ac.za
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Tue Sep 26 2017 - 06:00:02 PDT
