Just a comment. rdparm's "printLennardJones" command prints out the
pre-combined, A and B-coefficient arrays printed in the prmtop. If you
want to determine the original parameters as would be found in the frcmod
file, try the "printDetails" command of xparmed.py (GUI) or parmed.py
(ptraj-like text-based input).
All the best,
Jason
On Thu, Jul 5, 2012 at 5:48 PM, Ben Roberts <ben.roberts.geek.nz> wrote:
> 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
>
--
Jason M. Swails
Quantum Theory Project,
University of Florida
Ph.D. Candidate
352-392-4032
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu Jul 05 2012 - 15:30:02 PDT
