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
Received on Wed Jul 04 2012 - 17:30:03 PDT
