Hi Ben. Rdparm shows me that the only bond that exist with this
hydrogen is the one I have specified in the frcmod file (with the
hydroxyl oxygen). As to the electrostatics and van der Waals, I am not
exactly sure where would I obtain that information. Some help on this
would be great.
The apparent attempt to bond seems more like a very strong hydrogen
bond pulling the geometry towards that direction. It would not be what
one would expect for many reasons. Chemically speaking, the geometry
impediment would huge, since it is trying to from a three membered
ring. I've noticed, however, that the columns regarding this
particular dihedral angle (SG - ZN - O - Ho), when I ask in rdparm for
"dihedral :340.ZN" are like this:
Dihedral pk phase pn atoms
1166: 0.000 0.00 3.0 :31.SG :339.ZN :340.O :340.H1
(531,5311,5312,5313)
Is that normal?
Thank you
Fabrício Bracht
2012/7/4 Ben Roberts <ben.roberts.geek.nz>:
> Hi Fabrício,
>
> What does the program "rdparm" show in your topology? Is there any indication of a bond record or an angle record between the Ho and SG atoms? And are the electrostatics and van der Waals parameters appropriate, particularly when scaled to 1-4? And can you describe more of the connection between the torsional angle and the apparent attempt to "bond" - that is, are the bonds and angles involved chemically sensible?
>
> Cheers,
> Ben
>
> On 5/07/2012, at 6:20 AM, Fabrício Bracht wrote:
>
>> Hello. I have some interesting news. The procedure of addind the
>> glycam vdwradii seems to have solved the problem. I ran 10 ns using
>> pmemd.mpi and right after that I ran another 10 ns using pmemd.cuda.
>> Although the simulation now proceeds without any error, I have notice
>> that the dihedral angle for SG-ZN-OH-Ho (the cysteine Sulfur, the zinc
>> atom, the hydroxyl's oxygen and the hydroxyl's hydrogen) changes from
>> something around 130º-140º to -145º and than to 15º (more or less).
>> The visual aspect of the geometry is of a distorted one (for the
>> pmemd.cuda step). It looks as though the hydrogen is trying to bond
>> with the cysteine sulfur atom. Any help here would be great.
>> Thank you
>> Fabrício Bracht
>>
>> 2012/6/28 Fabrício Bracht <bracht.iq.ufrj.br>:
>>> Hello. Thank you for the advice. I've taken the GLYCAM path. Renamed
>>> my HO to Ho and started the process again. Will know if it worked in a
>>> few days. I'll contact the developers of GLYCAM in order to get more
>>> information, as well as some tips that my help my particular case.
>>>
>>> Thank you again for all the advice.
>>> Fabrício Bracht
>>>
>>> 2012/6/28 Ross Walker <ross.rosswalker.co.uk>:
>>>> Hi Fabricio,
>>>>
>>>> I suggest 'NOT' using parmed in this situation - if you start hacking the
>>>> prmtop file you will most likely forget what you did, or not interpret it
>>>> properly and when it comes time to publish and someone asks for the force
>>>> field parameters file and setup files to reproduce your work you won't have
>>>> them. Instead I would go back and modify the original parameterization files
>>>> you used in leap. This way it is clearly documented what you did and what
>>>> you changed. The hacking approach is useful to quickly test a hypothesis but
>>>> is no excuse for not being diligent in the approach in the first place.
>>>>
>>>> The VDW Radii do not come from the pdb file. They are VDW parameters
>>>> specified for the atom type in the parmXX.dat file or frcmod file. Here's an
>>>> example:
>>>>
>>>> http://ambermd.org/tutorials/advanced/tutorial1_adv/
>>>>
>>>> plc.frcmod:
>>>>
>>>> # modifications to force field for poplar plastocyanin
>>>>
>>>> MASS
>>>> SM 32.06
>>>> CU 65.36
>>>>
>>>> BOND
>>>> NB-CU 70.000 2.05000 #kludge by JRS
>>>> ...
>>>>
>>>> ANGLE
>>>> CU-NB-CV 50.000 126.700 #JRS estimate
>>>> ...
>>>>
>>>> DIHE
>>>> X -NB-CU-X 1 0.000 180.000 3.000
>>>> ...
>>>>
>>>> NONBON
>>>> CU 2.20 0.200
>>>> SM 2.00 0.200
>>>>
>>>> The values in the nonbond section correspond to the following:
>>>>
>>>> LTYNB , R , EDEP
>>>>
>>>> LTYNB Atom symbol.
>>>>
>>>> R The van der Waals radius of the atoms having the symbol
>>>> "LTYNB" (Angstoms)
>>>>
>>>> EDEP The 6-12 potential well depth. (kcal/mol)
>>>>
>>>> If you then look at $AMBERHOME/dat/leap/parm/parm99.dat you will see:
>>>>
>>>> H 0.6000 0.0157 !Ferguson base pair geom.
>>>> HO 0.0000 0.0000 OPLS Jorgensen,
>>>> JACS,110,(1988),1657
>>>> HS 0.6000 0.0157 W. Cornell CH3SH --> CH3OH FEP
>>>>
>>>> And there is your problem. The HO has no VDW interaction and thus can
>>>> collapse on other atoms. It doesn't if it is bonded or angled with them
>>>> since nonbond interactions are excluded for 1-2 and 1-3 terms but if it is a
>>>> 1-4 term or a 1-5 term or more than it has a charge interaction which could
>>>> be strong if they are different signs and pull the HO atom into it.
>>>>
>>>> One approach is what they do for hydroxyl's in Charmm All22:
>>>>
>>>> !V(Lennard-Jones) = Eps,i,j[(Rmin,i,j/ri,j)**12 - 2(Rmin,i,j/ri,j)**6]
>>>> !
>>>> !epsilon: kcal/mole, Eps,i,j = sqrt(eps,i * eps,j)
>>>> !Rmin/2: A, Rmin,i,j = Rmin/2,i + Rmin/2,j
>>>> !
>>>> !atom ignored epsilon Rmin/2 ignored eps,1-4 Rmin/2,1-4
>>>> ...
>>>> HT 0.000000 -0.046000 0.224500 ! ALLOW WAT
>>>>
>>>> Or you can look at what they do for GLYCAM_06H which is to define a Ho atom
>>>> type and give it a small VDW term:
>>>>
>>>> Ho 0.2000 0.0300 M.B. Tessier 2011, use by
>>>> permission only
>>>>
>>>> This is how they stop the hydroxyl hydrogens collapsing onto hydroxyl
>>>> oxygens in sugar rings.
>>>>
>>>> These are both hacks I believe, I don't think they actually fitted the
>>>> parameters here but you might want to ask the authors of those force fields
>>>> for details behind it. Then for your system you'll take your frcmod file
>>>> that you build during the parameterization and also the mol2 file (or lib
>>>> file) edit the atom types for the hydroxyl's to be something unused and then
>>>> provide the parameters (mass, bond, angle, dihedral) for this in the frcmod
>>>> file and finally at the end include the small VDW parameters that you want
>>>> for that atom.
>>>>
>>>> Good luck,
>>>>
>>>> All the best
>>>> Ross
>>>>
>>>>
>>>>> -----Original Message-----
>>>>> From: Fabrício Bracht [mailto:bracht.iq.ufrj.br]
>>>>> Sent: Thursday, June 28, 2012 1:29 PM
>>>>> To: AMBER Mailing List
>>>>> Subject: Re: [AMBER] NaN with GTX580 pmemd.cuda simulation with and
>>>>> without iwrap = 1
>>>>>
>>>>> Hello Ross. Thank you for the reply. I´ll restart the parameterization
>>>>> procedure. You suggest I create a new H type. I understand that
>>>>> setting a new type for this hydrogen would require me to include
>>>>> vdwradii for this atom. The problem is that I am not sure how I would
>>>>> do this. I would just have to rename the pdb source file? But then,
>>>>> how would I set the vdwradii? Can´t I just set the vdwradii for this
>>>>> particular atom in tleap, after having performed the parameterization?
>>>>> Wouldn´t the lmodprmtop program be of use in this case?.
>>>>> Thank you in advance
>>>>> Fabrício Bracht
>>>>>
>>>>> 2012/6/28 Ross Walker <ross.rosswalker.co.uk>:
>>>>>> Hi Fabricio,
>>>>>>
>>>>>> So you are running into one of the problems with the AMBER force
>>>>> field that
>>>>>> people have ignored for years and it was never too much of a problem
>>>>> since
>>>>>> people didn't typically run for long enough to see it. That is that
>>>>> hydroxyl
>>>>>> hydrogens have a zero VDW radii. In TIP3P water this is not a problem
>>>>>> because it is rigid and the H sits inside the VDW radii of the
>>>>> oxygen.
>>>>>> However, when you start moving to highly charged systems, a classic
>>>>> example
>>>>>> being phosphates you can end up with the H collapsing into one of the
>>>>> other
>>>>>> negatively charged atoms for which there is not a bond or angle term
>>>>> to
>>>>>> prevent it and as soon as that distance becomes very small so you
>>>>> divide by
>>>>>> zero and get NAN. The CPU is a little more tolerant of division by
>>>>> VERY
>>>>>> small numbers than the GPU and that's why it manifests itself more as
>>>>> a blow
>>>>>> up (trapped by vlimit) on the CPU and you get NAN on the GPU. Either
>>>>> way the
>>>>>> simulation is bust because this event is non physical. 1fs can help,
>>>>> since
>>>>>> the thing won't shoot off as much but it still won't really fix the
>>>>> problem.
>>>>>>
>>>>>> I would take a VERY careful look at the parameterization - that is
>>>>> the
>>>>>> underlying problem. To begin with you probably should not use type HO
>>>>> for
>>>>>> the hydrogen and instead create a new type. Then make sure this new H
>>>>> type
>>>>>> has at least a small VDW radii, a very small radii will help (look in
>>>>> the
>>>>>> frcmod file) but ideally one should probably attempt to refit the
>>>>> radii
>>>>>> properly for such cases but this is tricky. In your case I would just
>>>>> put a
>>>>>> small radii on the hydrogen (take a look at what the charmm people do
>>>>> with
>>>>>> TIP3P - which in charmm should really strictly be called TIP3P-CHARMM
>>>>> since
>>>>>> it isn't the original TIP3P parameterization due to the addition of
>>>>> small
>>>>>> radii on the hydroxyl hydrogens.)
>>>>>>
>>>>>> What this does highlight though is that it is important that people
>>>>> still
>>>>>> test things on the CPU and not just rely on the GPU code for
>>>>> everything. If
>>>>>> you see a problem with a GPU run then switch to CPU and analyse
>>>>> things
>>>>>> carefully. Don't expect the CPU code to also explode since it is more
>>>>>> tolerant of marginal stability cases but at least be aware that you
>>>>> might be
>>>>>> simulating something that is borderline stable at best.
>>>>>>
>>>>>> All the best
>>>>>> Ross
>>>>>>
>>>>>>> -----Original Message-----
>>>>>>> From: Fabrício Bracht [mailto:bracht.iq.ufrj.br]
>>>>>>> Sent: Thursday, June 28, 2012 12:09 PM
>>>>>>> To: AMBER Mailing List
>>>>>>> Subject: Re: [AMBER] NaN with GTX580 pmemd.cuda simulation with and
>>>>>>> without iwrap = 1
>>>>>>>
>>>>>>> Hello. Thank you for the reply. In fact, a time step of 1fs seems to
>>>>>>> work around the problem. At the same time as i was writing the last
>>>>>>> email, I had submitted a job with dt = 0.001 and it finished without
>>>>>>> any problems. I haven´t tested yet with pmemd.cuda, but I´ll do that
>>>>>>> later on.
>>>>>>> The two atoms involved in this case are O and H of the hydroxyl
>>>>>>> residue linked to the zinc atom of the active site. Where do I find
>>>>>>> the Lennard-Jones rmin values for them? The hole active site was
>>>>>>> parameterized following the steps on the MTK++ manual.
>>>>>>> I´ll follow your instructions as to getting a restart file from a
>>>>>>> previous step and test it out. I´ll inform you as soon as I have the
>>>>>>> results.
>>>>>>> Thank you
>>>>>>> Fabrício Bracht
>>>>>>>
>>>>>>> 2012/6/28 David A Case <case.biomaps.rutgers.edu>:
>>>>>>>> On Thu, Jun 28, 2012, Fabrício Bracht wrote:
>>>>>>>>
>>>>>>>>> Hi. I did what you asked. I ran the simulation again, in order to
>>>>>>> make
>>>>>>>>> sure that the error happened again. This time, it dies a lot
>>>>> earlier
>>>>>>>>> with the same warnings as last time (vlimit exceed....etc). I
>>>>>>>>> continued the simulation using as input the last restrt file
>>>>> written
>>>>>>>>> using the serial version of pmemd. The simulation crashed after a
>>>>>>> few
>>>>>>>>> hours , the last output lines are:
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> vlimit exceeded for step 112126; vmax = 91.1788
>>>>>>>>> vlimit exceeded for step 112176; vmax = 45.6806
>>>>>>>>>
>>>>>>>>> Coordinate resetting cannot be accomplished,
>>>>>>>>> deviation is too large
>>>>>>>>> iter_cnt, my_bond_idx, i and j are
>>>>>>> : 2 2690 5315 5316
>>>>>>>>>
>>>>>>>>
>>>>>>>> Thanks for the info. It looks to me like the problem has nothing
>>>>> to
>>>>>>> do (per
>>>>>>>> se) with GPU's, since you can reproduce it on a serial CPU using
>>>>>>> pmemd.
>>>>>>>> Unfortuantely, problems like this (generally with the force field)
>>>>>>> can be hard
>>>>>>>> to track down. Your mdin file looks fine to me, although it is
>>>>> safer
>>>>>>> to use
>>>>>>>> dt=0.001 rather than 0.002, especially for constant pressure
>>>>>>> simulations, and
>>>>>>>> if you system is not *quite* equilibrated.
>>>>>>>>
>>>>>>>> The information above indicates that the problem occurs around
>>>>> atoms
>>>>>>> 5315
>>>>>>>> and 5316. You might see what those are, as this might (or might
>>>>> not)
>>>>>>> give
>>>>>>>> a clue to the origin of the problem. Look hard (and/or let us
>>>>> know
>>>>>>> about)
>>>>>>>> any non-standard residues or molecules that are in your
>>>>>>> simulation. Do
>>>>>>>> you have any Lennard-Jones rmin values that are either zero or
>>>>> small
>>>>>>> (i.e.
>>>>>>>> less than typical values for that type of atom)? You may have to
>>>>> get
>>>>>>> a
>>>>>>>> restart file someplace a little before step 112126; (don't use
>>>>> ig=-1
>>>>>>> to do
>>>>>>>> this: copy the actual seed that is printed in your output file, so
>>>>>>> you will
>>>>>>>> get an identical trajectory.) Then you can set ntwx=1, and
>>>>> examine
>>>>>>> the
>>>>>>>> coordinates at every step as the system crashes....one hopes that
>>>>>>> this will
>>>>>>>> identify the source of the problem.
>>>>>>>>
>>>>>>>> ...good luck...dac
>>>>>>>>
>>>>>>>>
>>>>>>>> _______________________________________________
>>>>>>>> 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
>>>>>>
>>>>>>
>>>>>> _______________________________________________
>>>>>> 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
>>>>
>>>>
>>>> _______________________________________________
>>>> 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
>
>
> _______________________________________________
> 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
Received on Thu Jul 05 2012 - 14:00:02 PDT
