Hi Yun,
some dihedral torsions have a more complicated profile than a single cosine term alone could cover (else there is no point in using e.g. Ryckaert-Bellemans-potential). So In order to model a complicated dihedral torsion properly sometimes two or three terms with different dihedral periodicities are required. In that case the ff.dat shows entries with a negative sign at the periodicity value to prompt you/the parser to regard the next line as an additional contribution to the proper dihedral energy of this particular rotatable bond. The rdparm output apparently uses the "E" in front of the regular information for the same purpose. Consequently, the atoms involved appear repeatedly -- this does not imply the constants before are overridden, instead they are added to the ones in the line before.
I could not reproduce your troubles in creating top/crd-files with a modified force field version. With your file (mono.pdb) as input I can create prmtop + inpcrd files for the complete compound with hydrogens. Using GLYCAM_06c.dat, this gives me a topology file equivalent to mono.prmtop (except for the OS-CG-OS-CG torsions, where values apparently have been updated compared to the version you used ). I then made a copy of GLYCAM_06c.dat (attached as a reference) and changed the parameters for OH-CG-CG-OH according to our previous thoughts (glycam_phase_shift.dat, attached) and can perfectly well create a top-file (with the leapin file attached). Consequently I can create two variants of gromacs-topology files which differ in the RB-coefficients for one single dihedral type (occurring two times) with a plenty error messages, but nothing fatal. The resulting files appear formally correct.
To illustrate the consequences of using negative force constants in the input-ff.dat I attached the diagram dihedral_profiles.eps, created with gpt_in (look there for the equations and force constants used). Admittedly the effect of the bug in amb2gmx.pl is not too dramatic in this case -- but it was a small force constant of 0.1 kcal in this example, for other torsions this issue may have more impact.
Hope that helps
Astrid
----- Ursprüngliche Mail -----
Von: "Yun Shi" <yunshi09.gmail.com>
An: "AMBER Mailing List" <amber.ambermd.org>
Gesendet: Samstag, 10. September 2011 01:51:48
Betreff: Re: [AMBER] how to fill in missing parameters in tleap?
Hi Astrid,
I changed the code, and converted a simple monosaccharide. But I noticed
that for eg:
............
73: -1.100 0.00 1.0 :2.O5 :2.C5 :2.C4 :2.O4 (8,9,15,17)
E 74: 0.250 0.00 2.0 :2.O5 :2.C5 :2.C4 :2.O4 (8,9,15,17)
...............
So for the same sequence (atoms 8, 9, 15, 17), there are two set of dihedral
parameters, one of which is negative, which I am afraid would be
problematic. What that "E" in front of "74" mean?
Anyways, I am having trouble with tleap now after playing with .dat file,
something like
saveamberparm mono monom.prmtop monom.inpcrd
Checking Unit.
Building topology.
Building atom parameters.
For atom: .R<OME 1>.A<H1 1> Could not find type: H1
For atom: .R<OME 1>.A<CH3 2> Could not find type: CG
For atom: .R<OME 1>.A<H2 3> Could not find type: H1
For atom: .R<OME 1>.A<H3 4> Could not find type: H1
For atom: .R<OME 1>.A<O 5> Could not find type: OS
For atom: .R<0hA 2>.A<C1 1> Could not find type: CG
For atom: .R<0hA 2>.A<H1 2> Could not find type: H2
For atom: .R<0hA 2>.A<O5 3> Could not find type: OS
For atom: .R<0hA 2>.A<C5 4> Could not find type: CG
For atom: .R<0hA 2>.A<H5 5> Could not find type: H1
For atom: .R<0hA 2>.A<C6 6> Could not find type: CG
For atom: .R<0hA 2>.A<H61 7> Could not find type: H1
For atom: .R<0hA 2>.A<H62 8> Could not find type: H1
For atom: .R<0hA 2>.A<H63 9> Could not find type: H1
For atom: .R<0hA 2>.A<C4 10> Could not find type: CG
For atom: .R<0hA 2>.A<H4 11> Could not find type: H1
For atom: .R<0hA 2>.A<O4 12> Could not find type: OH
For atom: .R<0hA 2>.A<H4O 13> Could not find type: HO
For atom: .R<0hA 2>.A<C3 14> Could not find type: CG
For atom: .R<0hA 2>.A<H3 15> Could not find type: H1
For atom: .R<0hA 2>.A<O3 16> Could not find type: OH
For atom: .R<0hA 2>.A<H3O 17> Could not find type: HO
For atom: .R<0hA 2>.A<C2 18> Could not find type: CG
For atom: .R<0hA 2>.A<H2 19> Could not find type: H1
For atom: .R<0hA 2>.A<O2 20> Could not find type: OH
For atom: .R<0hA 2>.A<H2O 21> Could not find type: HO
Parameter file was not saved.
showed up when I tried to save top files, and even though I changed the .dat
file back to original, it's still not working.
Attached is the .pdb, .prmtop, .inpcrd, and top and gro file generated by
MODIFIED amb2gmx.pl. Maybe someone could try the approach you mentioned
earlier before I got my tleap working.
Yun
On Fri, Sep 9, 2011 at 4:24 AM, Astrid Maaß <astrid.maass.scai.fraunhofer.de
> wrote:
> Hi Yun,
>
> you have outlined the procedure the way I would do it.
>
> You already have identified an "offensive" set of atoms with rdparm. In
> your example these atoms are
> 3.H2N :3.N2 :3.C2N :3.O2N. The residue number (":3" in this case) is
> not relevant, but the atom names follow the leading '.'-character are the
> ones you need. In glycam04.prep I find these names e.g. for the residue 1YB,
> where the atom names and types (plus their connectivity) are listed. So you
> need the atom-type combination "H -N -C -O" (or possibly the reversed
> pattern; mind the blanks). In fact, for this torsion, glycam04.dat shows the
> following entry (where, ironically, you want to reverse the correction ;-)
> ).
>
> H -N -C -O 1 2.00 0.0 -1. Parm94 -
> Corrected to remove phase shift
> 1 -2.50 0.0 2.
>
> Modify as suggested earlier and proceed as before.
>
> Good luck!
>
> Astrid
>
>
>
> ----- Ursprüngliche Mail -----
> Von: "Yun Shi" <yunshi09.gmail.com>
> An: "AMBER Mailing List" <amber.ambermd.org>
> Gesendet: Freitag, 9. September 2011 08:47:16
> Betreff: Re: [AMBER] how to fill in missing parameters in tleap?
>
> Hi Astrid,
>
> I would be very happy to share.
>
> But before I try this, I just want to make sure if I could just find the
> pattern of atomtypes concerned (with negative dihedral force constants as
> displayed by rdparm) in Glycam_06g.dat file, modify corresponding
> parameters
> in the way you mentioned before, and then use tleap/xleap to generate
> prmtop
> again, which is subsequently converted by amb2gmx.pl?
>
> Thanks,
>
> Yun
>
> On Thu, Sep 8, 2011 at 11:24 PM, Astrid Maaß <
> astrid.maass.scai.fraunhofer.de> wrote:
>
> > Hi Yun,
> >
> > ----- Ursprüngliche Mail -----
> > > Von: "Yun Shi" <yunshi09.gmail.com>
> > > An: "AMBER Mailing List" <amber.ambermd.org>
> > > Gesendet: Donnerstag, 8. September 2011 22:50:55
> > > Betreff: Re: [AMBER] how to fill in missing parameters in tleap?
> > >
> > > Hi Astrid,
> > >
> > > Thanks for the advice. The problem in amb2gmx.pl has not been fixed,
> > > which
> > > still do not recognize negative pK. So what you were saying is to
> > > change
> > > something like:
> > >
> > > E 53: -2.500 0.00 2.0 :3.H2N :3.N2 :3.C2N :3.O2N
> > > (45,44,46,47)
> > >
> > > to
> > >
> > > E 53: 2.500 180.00 2.0 :3.H2N :3.N2 :3.C2N :3.O2N
> > > (45,44,46,47)
> > >
> > > ?
> > >
> > > And then use amb2gmx.pl to convert?
> > >
> >
> > yes, this is exactly, what I meant.
> >
> >
> > > I also wonder if it's OK to change the ">" to "!=" as in
> > > ...........................
> > >
> > > # get all force constants for each line of a dihedral #
> > > my $lines = $i -1 +$numijkl;
> > > for(my $j=$i;$j<=$lines;$j++){
> > > my $period = abs($pn{$j});
> > > if($pk{$j}>0) {
> > > $V[$period] = 2*$pk{$j}*$cal/$idivf{$j};
> > > }
> > >
> > > ...........................
> > >
> > > Regards,
> > >
> > > Yun
> > >
> >
> > Adressing the problem by fixing the code is certainly better, but since I
> > am not much of a programming person, I am not sure, whether the fix you
> > suggest below remedies the problem. If you decide to try both variants
> > (modified input FF + old code; old FF + modified code), I think the
> > resulting gromacs topologies should finally display consistent
> > RB-coefficients.
> >
> > I would appreciate if you could share the outcome of this with us.
> >
> > Regards,
> >
> > Astrid
> >
> >
> > >
> > > On Wed, Sep 7, 2011 at 12:14 AM, Astrid Maaß <
> > > astrid.maass.scai.fraunhofer.de> wrote:
> > >
> > > > Hi Yun,
> > > >
> > > >
> > > > as per amber file format specifications, the minus indicates that
> > > > the
> > > > dihedral potential of this particular set of atoms is composed of
> > > > more than
> > > > one cosine term; this appears to be properly recognized, I had not
> > > > troubles
> > > > with amb2gmx.pl in this context.
> > > >
> > > > The numbers in columns 6 - 11 of the gromacs-topology are the
> > > > Ryckaert-Bellemans-constants you get AFTER converting the
> > > > amber-topology to
> > > > the gromacs. Do not worry about the signs in the GROMACS-topology.
> > > >
> > > > What you should check with rdparm, if you create a topology file
> > > > for input
> > > > in amb2gmx.pl with tleap/xleap/sleap anyways, whether dihedrals
> > > > with
> > > > negative force constants are present in the AMBER-topology file.
> > > > The rdparm
> > > > output is easier to read than the topology file itself and tells
> > > > you, which
> > > > torsion exactly i.e. which atoms are involved and thus which
> > > > pattern of
> > > > atomtypes is concerned (you can check the link between an atom's
> > > > name and
> > > > type in prep or mol2 files).
> > > >
> > > > If you find suspect torsions/force constants in the rdparm-output
> > > > (or the
> > > > FF-input.dat itself), create a copy of the FF-input.dat and modify
> > > > within
> > > > this copy the respective force constants & phase shifts for the
> > > > identified
> > > > patterns. Then create with the modified FF-input.dat a new set of
> > > > top- and
> > > > crd-files for conversion with amb2gmx.pl. As I said before, I am
> > > > not sure
> > > > whether the problem might have been fixed since the time when I
> > > > encountered
> > > > it, so check, whether the whole fuss is necessary and if the old
> > > > and the new
> > > > gromacs-topologies are indeed different.
> > > >
> > > > Good luck!
> > > >
> > > > Astrid
> > > >
> > > >
> > > >
> > > >
> > > >
> > > >
> > > >
> > > > ----- Ursprüngliche Mail -----
> > > > Von: "Yun Shi" <yunshi09.gmail.com>
> > > > An: "AMBER Mailing List" <amber.ambermd.org>
> > > > Gesendet: Dienstag, 6. September 2011 20:22:06
> > > > Betreff: Re: [AMBER] how to fill in missing parameters in tleap?
> > > >
> > > > Hi Astrid,
> > > >
> > > > I saw something like:
> > > >
> > > > [ dihedrals ]
> > > > ;i j k l func C0 ... C5
> > > > 19 18 20 21 3 0.75312 2.25936 0.00000
> > > > -3.01248 0.00000 0.00000 ;
> > > > 41 1 18 19 3 0.20920 0.62760 0.00000
> > > > -0.83680 0.00000 0.00000 ;
> > > >
> > > > after converting glycam ff-parametized saccharide to gromacs
> > > > topology
> > > > format. I referred to gromacs manual, and knew that '3' in the
> > > > fifth column
> > > > indicates the Ryckaert-Bellemans-potential you mentioned. But could
> > > > you
> > > > tell
> > > > me what those numbers in columns 6 to 11 mean? Which is force
> > > > constant and
> > > > which is phase shift?
> > > >
> > > > Thanks,
> > > >
> > > > Yun
> > > >
> > > >
> > > > On Tue, Sep 6, 2011 at 8:21 AM, Yun Shi <yunshi09.gmail.com> wrote:
> > > >
> > > > > Hi Astrid,
> > > > >
> > > > > Thank you for the info.
> > > > >
> > > > > But in addition to negative force constant, there are also some
> > > > > negative
> > > > > multiplicity for GLYCAM dihedral parameters. So would these
> > > > > negative
> > > > > multiplicities be a problem for amb2gmx as well?
> > > > >
> > > > > Thanks,
> > > > >
> > > > > Yun
> > > > >
> > > > >
> > > > >
> > > > > On Tue, Sep 6, 2011 at 2:07 AM, Astrid Maaß <
> > > > > astrid.maass.scai.fraunhofer.de> wrote:
> > > > >
> > > > >> Hi Yun Shi,
> > > > >>
> > > > >>
> > > > >> I used amb2gmx.pl some time ago and if I remember correctly, for
> > > > >> torsions, amb2gmx.pl converts the Amber-cosine series in a
> > > > >> Ryckaert-Bellemans-potential. This was OK in most of the cases,
> > > > >> but when
> > > > >> amb2gmx.pl encounters a negative force constant it may not
> > > > >> translate
> > > > the
> > > > >> potential correctly. I do not know, whether this behavior has
> > > > >> been fixed
> > > > in
> > > > >> the meantime, but you could bypass it by using a phase shift of
> > > > >> 180
> > > > instead
> > > > >> of a negative sign in the force constant for the respective
> > > > >> term. You
> > > > might
> > > > >> use rdparm to check whether this difficulty would apply to your
> > > > >> system
> > > > >> anyways.
> > > > >>
> > > > >> Hope that helps
> > > > >> Astrid
> > > > >>
> > > > >>
> > > > >>
> > > > >>
> > > > >> ----- Ursprüngliche Mail -----
> > > > >> Von: "Yun Shi" <yunshi09.gmail.com>
> > > > >> An: "AMBER Mailing List" <amber.ambermd.org>
> > > > >> Gesendet: Dienstag, 6. September 2011 01:22:19
> > > > >> Betreff: Re: [AMBER] how to fill in missing parameters in tleap?
> > > > >>
> > > > >> Hi all,
> > > > >>
> > > > >> For this specific molecule, I wonder if it's possible to
> > > > >> manually set
> > > > >> parameters for the non-standard parts.
> > > > >>
> > > > >> I was thinking about using GLYCAM and amber99SB, while avoiding
> > > > >> the use
> > > > of
> > > > >> GAFF force field, with the help of amb2gmx.pl and GROMACS
> > > > >> software.
> > > > >>
> > > > >> 1. Set geometrical parameters and atomic charges to all standard
> > > > residues
> > > > >> using GLYCAM and 99SB.
> > > > >>
> > > > >> 2. Derive RESP charges for the modified TRP and thio-glycosidic
> > > > >> linkage.
> > > > >>
> > > > >> 3. With GROMACS topology files, set geometrical parameters for
> > > > >> the
> > > > >> modified
> > > > >> TRP and thio-glycosidic linkage according to 99SB. That is,
> > > > >> define atom
> > > > >> types, bonds (CT-S), angles (CT-S-CT, CT-CT-S), and dihedrals
> > > > >> (X-CT-S-X)
> > > > >> with corresponding parameters from 99SB. And for N, H, CA, CB,
> > > > >> and
> > > > indole
> > > > >> parts of the modified TRP, just use whatever geometrical
> > > > >> parameters from
> > > > >> the
> > > > >> standard TRP in 99SB.
> > > > >>
> > > > >> 4. Combine all these parts in GROMACS format.
> > > > >>
> > > > >>
> > > > >> On Sun, Aug 28, 2011 at 6:12 AM, case <case.biomaps.rutgers.edu>
> > > > >> wrote:
> > > > >>
> > > > >> > On Fri, Aug 26, 2011, Jason Swails wrote:
> > > > >> > >
> > > > >> > > Use parmchk. If you run the program with no arguments, it
> > > > >> > > tells you
> > > > >> how
> > > > >> > to
> > > > >> > > use it. By default, it will search gaff.dat and pull out
> > > > >> > > reasonable
> > > > >> > > parameters, dumping them to an frcmod file. If there's
> > > > >> > > anything it
> > > > >> > simply
> > > > >> > > *can't* find, it will put a placeholder (all parameters are
> > > > >> > > 0) and
> > > > say
> > > > >> > > "ATTN, need revision". You have to fix those parameters.
> > > > >> >
> > > > >> > Parmchk does a bit more than just look in gaff.dat: it
> > > > >> > estimates
> > > > missing
> > > > >> > parameters by several algorithms. In its default mode,
> > > > >> > parameters in
> > > > >> > gaff.dat
> > > > >> > are not put into the frcmod file; parameters it estimates on
> > > > >> > its own
> > > > are
> > > > >> > output; finally, things it can't figure out get placeholders.
> > > > >> > >
> > > > >> > > However, I would warn you that gaff dihedrals should not
> > > > >> > > necessarily
> > > > >> be
> > > > >> > > accepted. Run short simulations to make sure they make
> > > > >> > > sense, as
> > > > you
> > > > >> may
> > > > >> > > have to modify some of them by hand.
> > > > >> >
> > > > >> > This is excellent advice.
> > > > >> >
> > > > >> > ...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
>
_______________________________________________
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
- application/octet-stream attachment: gpt_in
- application/octet-stream attachment: leapin
Received on Mon Sep 12 2011 - 06:30:04 PDT
