Hi Fabrício,
The electrostatics and van der Waals parameters can also be obtained using rdparm. There's a command one can use (called, I think, "help") that prints out a list of commands one can use to get a list of the electrostatic and Lennard-Jones parameters.
The torsional parameter you provided removes any remaining doubt in my mind that torsions are *not* directly responsible for what you're seeing. I had already suspected that was the case. A torsional parameter, being solely concerned with rotation around the central bond, should not pull the atoms in 1- and 4- positions towards each other, as that would involve a force component orthogonal to any "torsional" force.
To the best of my knowledge, attractive forces between two particles, of the sort that would cause unusual geometries, could only result from:
- A spurious bond record;
- A harmonic restraint; or
- Non-bonding (i.e., electrostatic and/or van der Waals) interactions.
You've already said there is only one bond to that hydrogen, which rules out the spurious bond explanation. Likewise, I would guess you aren't using any harmonic restraints, though it would be helpful if you could confirm that. Which leaves only the non-bonding interactions, that I've already alluded to.
Thus, it seems likely to me that you have highly charged particles and the combination of van der Waals parameters and bond/angle/torsion parameters isn't enough to keep them apart. Given that, perhaps you could supply for us more bits of information:
- The atomic partial charges of the hydrogen and sulphur atoms
- The Lennard-Jones parameters of those atoms
- The bond parameters of SG-ZN, ZN-OH and OH-Ho
- The angle parameters of SG-ZN-OH and ZN-OH-Ho
- A list of which particular bonds and angles among these are significantly distorted from their reference values (r0 and theta0) as the hydrogen and sulphur atoms approach each other
- A description of which parameter set you used. Did you make your parameters using MCPB? If so, did you use the MD (Z-matrix) method or the seminario method to get your parameters? I ask because one of these creates overly weak bonds, and the other overly weak angles, so I've found it necessary to mix the two in the past.
Cheers,
Ben
On 6/07/2012, at 8:42 AM, Fabrício Bracht wrote:
> 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
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu Jul 05 2012 - 15:00:02 PDT
