Re: [AMBER] Using the RED Server Dev. to parameterize a model transition state

From: FyD <fyd.q4md-forcefieldtools.org>
Date: Fri, 13 Dec 2013 16:30:51 +0100

Dear Kamali,

> I also agree that I shouldn't use these results for my final calculations;
> I just used the OPT_Calc = "OFF1" option as part of a "practice run" so
> that I could see how to successfully submit a job to R.E.D Python and how
> the outputs would look.
>
> In a previous message you also asked why I used the INTRA-MCC constraints
> that I did: this was also part of the "practice". I chose atoms to
> constrain that made sense to me, and I just wanted to see what would
> happen.

You constrain atoms (using intra-mcc) that needs to be; for instance
to generate a molecular fragment from a whole molecule; in your case
you study a TS so I do not think you need to generate a molecular
fragment (except if you do want to connect it to a larger entity).

> I actually don't think I need to constrain any atoms beyond the
> default setting; is it possible to also use defaults for INTRA-MCC?

No

> Maybe I should start a different thread for this question, but I have
> actually tried to submit job P8902 with OPT_Calc = OFF and providing my own
> QM output. However, I've gotten errors in the R.E.D.-
> Server-14651.master0.q4md-forcefieldtools.org file:

ok I looked at the minimum you provided in the P8902 job. This is not
the TS you are interested in... please compare your PDB input file and
this QM output i.e. the last set of Cartesian coordinates...

> *...*
> *Reading JOB1*.log for molecule : 1*
> *Find "Normal termination of Gaussian"*
> *Find "Normal termination of Gaussian"*
> *Find "Stationary point found"*
> *Find "Standard orientation"*
> *Element in JOB1 is not the same with PDB*
> *Find big error in log : /home/u0263B7BZ/P8902/JOB1-gau_m1-1.log*
> *JOB1*.log is not correct , can not continue......*
> *Information of end*
> *Time : Wed Dec 11 22:26:46 CET 2013*
>
> But to me the JOB1-gau_m1-1.log looks fine. Did I have to modify my QM
> output (similar to what one must do for R.E.D Server Perl)?

try removing the Frequency job from the Mol_red1.log file as described
in red color at:
http://q4md-forcefieldtools.org/Tutorial/Tutorial-1.php#6
(see -4-The last step consists of checking and preparing the P2N and
QM output files for the R.E.D. III.x execution)

>> In these conditions _the key_ is the definition of the atom
>> connectivities from this PDB file (atom typing & FF parameters are
>> based on these atom connectivities). You have C-H bonds with 1.5 Angs
>> length in your PDB input file. Are you sure an atom connectivity is
>> created in this case (VMD does not...)?
>
> You're correct; those bond lengths should be fixed, and I will do that
> before the next time I submit my structure.
>
>
>> Your best bet here (because RED Python computes atom connectivities
>> after the geometry optimization step and you skip this step...) is to
>> provide the atom connectivities *YOU* want in the PDB input file and
>> *prevent* their computation.
>>
>> See
>> http://q4md-forcefieldtools.org/REDS-Development/Demo2-Files/Project.config
>> MOLECULE1-CALCONECT = OFF
>>
>> then you need to carefully think which atom types you want: obviously
>> two CT-H1 or CT-H2 bonds (CT, H1, H2 are atom types) cannot have two
>> different equilibrium values; i.e. ~1.0 and ~1.5 Angs -> a new atom
>> type has to be created...
>>
>> Moreover your phosphor atom has 5 bonds I would set up a new atom type
>> for it...
>
> Thank you for the suggestion, Francois. But I am wondering if the new
> phosphorus atom type is necessary? I ask because it is possible in nature
> for phosphorus to form 5 bonds.

look at the atom type definition in the Amber FF for the phosphorus
atom in the phosphate group; it is related to this NA structure:
http://en.wikipedia.org/wiki/File:DNA_chemical_structure.svg

>> Finally I really hope the Cartesian coordinates you have were checked
>> by QM: once again better providing a QM output (with one imaginary
>> frequency) as input for RED Python...
>
> For job P8902, I used Gaussian09 with these options:
>
> #P hf/6-31G* Opt=(Tight,CalcFC) Freq SCF(Conver=8) Test
>
> I tried to answer this question on my own, but I got confused: why do I
> need imaginary frequencies for my calculations, and why only just one?

See http://www.gaussian.com/g_tech/g_ur/k_opt.htm
  & use the QST2 or QST3 keyword.

one imaginary frequency for a TS is what I learned: I think you need
to read about that ;-)

regards, Francois


>> Quoting FyD <fyd.q4md-forcefieldtools.org>:
>>
>> > Kamali,
>> >
>> > please check the PDB file format you use as input...
>> > see
>> http://q4md-forcefieldtools.org/REDS-Development/Demo2-Files/readme.txt
>> >
>> > regards, Francois
>> >
>> >
>> > Quoting FyD <fyd.q4md-forcefieldtools.org>:
>> >
>> >> Dear Kamali,
>> >>
>> >>> I have been trying a few things out with my transition state structure
>> of
>> >>> RNA cleavage, and I had a few follow-up questions. My questions are
>> >>> regarding my job P8824:
>> >>
>> >> ok
>> >>
>> >> in the PDB input file:
>> >> ATOM 15 O2P RT0 1 31.511 16.369 54.925
>> >> O
>> >> ATOM 16 O5' RT0 6 32.667 14.580 56.630
>> >> O
>> >> ATOM 17 P RT0 1 30.934 15.202 55.888
>> >> P
>> >> ->
>> >> ATOM 15 O2P RT0 1 31.511 16.369 54.925
>> >> O
>> >> ATOM 16 O5' RT0 1 32.667 14.580 56.630
>> >> O
>> >> ATOM 17 P RT0 1 30.934 15.202 55.888
>> >> P
>> >>
>> >> & look at the residue name generated by R.E.D. Python
>> >>
>> >> You deactivated many keywords in the Project.config file letting the
>> >> defaults to be executed; this is correct. I am glad to see you can do
>> >> the job without a tutorial. We need to provide a tutorial. The
>> >> description of the keywords in the Project.config and Configuration.py
>> >> files is too limited...
>> >>
>> >>> Following your suggestion in the private communication (when we first
>> >>> talked about this structure), for this job I just ran the PDB through
>> the
>> >>> RED Server Dev. with the option OPT_Calc = Off1 in the Project.config
>> file,
>> >>> as a practice run.
>> >>
>> >> OPT_Calc = Off1 means you skip the QM job and use the Cartesian
>> >> coordinates from the PDB file to be used in MEP computation. (I am
>> >> really not a big fan of this approach; I would prefer you provide a QM
>> >> job for a Cartesian coordinate set with one imaginary frequency). ok
>> >> for a 'dirt cheap' job ;-) (just my personal opinion)
>> >>
>> >>> I was wondering if there is any documentation for the
>> >>> files generated by RED Python?
>> >>
>> >> Not yet - so far we work on a new version of R.E.D. Python with new
>> >> features...
>> >>
>> >> Just ask for help; sorry for that...
>> >>
>> >>> I found two types of files in the
>> >>> P8824/Data-R.E.D.Server/Mol_m1 folder: *sm* files (which I know are
>> files
>> >>> containing intra-molecular charge constraints) and the *ia*: what does
>> the
>> >>> *ia* mean? I also have two .mol2 files: Mol-ia0_m1-c1.mol2
>> >>
>> >> sm means files related to single molecule charge fit (without intra-mcc)
>> >> ia means files related to charge fitting carried out with intra-mcc
>> >> (because you use 'MOLECULE1-INTRA-MCC1 = 0.0 | 1 4 11 | Keep')
>> >> mm means files related to multiple molecule fit (i.e. number of
>> >> molecules > 1)
>> >>
>> >>> and Mol-sm_m1-c1.mol2, but these files appear to have two columns of
>> >>> charges rather than one next to the rightmost **** column. Could you
>> tell
>> >>> me which column is the correct one for charges, or alternatively point
>> me
>> >>> to documentation that answers my questions?
>> >>
>> >> See
>> >>
>> http://cluster.q4md-forcefieldtools.org/~x/P8824/Data-R.E.D.Server/Mol_m1/Mol-sm_m1-c1.mol2
>> >>
>> >> 1 C4' 30.706000 13.015000 52.944000 CT 1 U01 0.0520 0.8780
>> ****
>> >>
>> >> 0.0520 is the column of the atomic charge value
>> >> 0.8780 is the column of the atomic polarizability
>> >> (useful for the non-additive FF model)
>> >>
>> >> ---
>> >>
>> >> I do not understand why you use: MOLECULE1-INTRA-MCC1 = 0.0 | 1 4 11 |
>> Keep
>> >> in
>> >>
>> http://cluster.q4md-forcefieldtools.org/~x/P8824/Data-R.E.D.Server/Mol_m1/Mol-ia0_m1-c1.mol2
>> >>
>> >> i.e. 0.6957-0.5092-0.1865=0.0000
>> >>
>> >> then look at the impact of the INTRA-MCC1:
>> >>
>> http://cluster.q4md-forcefieldtools.org/~x/Project/P8824/Data-R.E.D.Server/Mol_m1/Statistics_m1.txt
>> >>
>> >> See the table:
>> >> Charge values for single molecule fit; without (SM) vs with (IA)
>> INTRA-MCC:
>> >> -> the impact seems quite important...
>> >>
>> >> You could also look at :
>> >> - the RRMS values of the corresponding fits.
>> >> - the es*.pdb files to display where the errors are.
>> >> See http://q4md-forcefieldtools.org/RED/resp/#intro
>> >> espdb
>> >> -j PDB-like file with relative residual in the TempFactor field
>> >> esqpotpdb
>> >> -y PDB-like file with input MEP values in the TempFactor field
>> >> (atomic units)
>> >> esmpotpdb
>> >> -z PDB-like file with MEP values for new charges in TempFactor field
>> >> (atomic units)
>> >>
>> >> let us know if you need more information...
>> >>
>> >> regards, Francois
>> >>
>> >>
>> >>> On Tue, Dec 10, 2013 at 8:12 AM, Kamali Sripathi
>> >>> <ksripath.umich.edu> wrote:
>> >>>
>> >>>> Dear Francois,
>> >>>>
>> >>>> Thank you very much for the detailed information, and I'm sorry for
>> the
>> >>>> confusion. I would have mentioned that I had questions about a topic
>> >>>> related to our private discussion, but I didn't want you to have
>> >>>> to look up
>> >>>> that communication, and so I tried to start a new thread :)
>> >>>> Looking back, I
>> >>>> see that my initial email in this thread was unfortunately less clear
>> than
>> >>>> I intended.
>> >>>>
>> >>>> I will try the protocol you suggested, but I think it would be better
>> for
>> >>>> me to use the RESP-A1 charge type because I'm working with RNA rather
>> than
>> >>>> carbohydrates, and using AMBER10 for my simulations.
>> >>>>
>> >>>> Thank you very much again, and have a great day,
>> >>>>
>> >>>> Kamali
>> >>>>
>> >>>>
>> >>>> On Tue, Dec 10, 2013 at 2:52 AM, FyD <fyd.q4md-forcefieldtools.org>
>> wrote:
>> >>>>
>> >>>>> Dear Kamali,
>> >>>>>
>> >>>>> > No, thank you very much for the offer. I have the PDF and refer to
>> it
>> >>>>> often
>> >>>>> > because there are concepts that I'm still learning :)
>> >>>>>
>> >>>>> Yes a looot to learn in this work...
>> >>>>>
>> >>>>> > I have two final question: from examples on the RED Server (perhaps
>> >>>>> this is
>> >>>>> > not applicable to RED Perl or RED Python), it appears that one can
>> also
>> >>>>> set
>> >>>>> > the INTRA-MCC constraints to defined values, e.g., sugar carbons to
>> >>>>> > standard values in the AMBER force field instead of setting the
>> sum of
>> >>>>> > certain atoms to 0 or to another value. Do you suggest one way
>> over the
>> >>>>> > other?
>> >>>>>
>> >>>>> I would constrain hydrogen atoms (only these in the methylene &
>> methyl
>> >>>>> groups, using the 'F' flag in the intra-mcc to keep them in the FF
>> >>>>> lib) _only_ if you use the GLYCAM force field. So the 'RESP-C2'
>> charge
>> >>>>> model (to be provided in the Configuration.py file to overwrite the
>> >>>>> default RESP-A1):
>> >>>>> see
>> http://q4md-forcefieldtools.org/REDS-Development/popup/popkeyword.php
>> >>>>>
>> >>>>> Here once again I would use INTRA-MCC1; in general one always uses
>> >>>>> INTRA-MCC1 except when one realizes that some atoms are not
>> >>>>> equivalenced as they should be ;-)
>> >>>>>
>> >>>>> > I understand that it very much depends on one's case, but I was
>> >>>>> > wondering if you had any thoughts. I was thinking that I might
>> >>>>> constrain
>> >>>>> > the charges of some of the sugar atoms that are 5' of the
>> >>>>> transition-state
>> >>>>> > pentacovalent phosphate to their values in the AMBER force fields
>> >>>>> instead
>> >>>>> > of summing their charges to 0.
>> >>>>>
>> >>>>> oh oh I got it: You continue the discussion from the private
>> >>>>> assistance...
>> >>>>> Not easy to follow you. In your case I would:
>> >>>>> -1 characterize the TS using the QM program of your choice
>> >>>>> -2 run RED Python using CHR_TYP = "RESP-C2", OPT_Calc = "Off" &
>> >>>>> MEPCHR_Calc = "On" using the QM log file & PDB file as inputs
>> >>>>> -3 study what you get: here you will be able to:
>> >>>>> - check the atom connectivities generated in the mol3 file
>> >>>>> - re-run RED Python using Re_Fit = "On" (i.e. upload/provide the
>> >>>>> entire/previous RED Python job in the archive)
>> >>>>> - provide new FF parameters in the frcmod.user (see the
>> >>>>> frcmod.unknown file generated in the previous job) in the archive
>> file
>> >>>>> See
>> >>>>>
>> http://q4md-forcefieldtools.org/REDS-Development/RED-Server-demo1.php
>> >>>>>
>> >>>>>
>> http://q4md-forcefieldtools.org/REDS-Development/popup/poptestcase.php
>> >>>>>
>> >>>>>
>> http://q4md-forcefieldtools.org/REDS-Development/Demo2-Files/frcmod.user
>> >>>>> - provide new FF atom types in the Project.config file for key
>> atom
>> >>>>> center(s)
>> >>>>> MOLECULE1-ATMTYPE = ...
>> >>>>> (here you can load your own FF as input to R.E.D. Python)
>> >>>>>
>> >>>>> When ready just request a private assistance and provide the PXXXX
>> >>>>> R.E.D. Server job name in the body of your email.
>> >>>>>
>> >>>>> > I lastly wanted to double-check that the INTER-MCC and INTER-MEQ
>> >>>>> keywords
>> >>>>> > are only necessary if one is deriving charges for more than one
>> >>>>> molecule,
>> >>>>> > and so would not be applicable in my one-omolecule one-conformation
>> >>>>> case?
>> >>>>>
>> >>>>> Yes INTRA means within a molecule; INTER means between molecules.
>> >>>>>
>> >>>>> regards, Francois
>> >>>>>
>> >>>>>
>> >>>>> >> regards, Francois
>> >>>>> >>
>> >>>>> >>
>> >>>>> >> > On Sat, Dec 7, 2013 at 4:30 PM, FyD <
>> fyd.q4md-forcefieldtools.org>
>> >>>>> >> wrote:
>> >>>>> >> >
>> >>>>> >> >> Dear Kamali,
>> >>>>> >> >>
>> >>>>> >> >> > I am using the RED Server Development to parameterize a model
>> >>>>> compound
>> >>>>> >> >> for
>> >>>>> >> >> > the transition state and different protonation states of one
>> or
>> >>>>> both
>> >>>>> >> of
>> >>>>> >> >> the
>> >>>>> >> >> > nonbridging oxygens in ribozyme cleavage.
>> >>>>> >> >>
>> >>>>> >> >> Classical geometry optimization as implemented in R.E.D.
>> >>>>> will likely
>> >>>>> >> >> generate a true minimum; not a TS.
>> >>>>> >> >>
>> >>>>> >> >> So better first characterizing the TS & loading the
>> >>>>> corresponding QM
>> >>>>> >> >> log file as input by setting OPT_Calc = "Off" & MEPCHR_Calc =
>> "On"
>> >>>>> in
>> >>>>> >> >> the Configuration.py file
>> >>>>> >> >>
>> >>>>> >> >> > My questions are regarding a
>> >>>>> >> >> > portion of the Project.config file that one may optionally
>> >>>>> include to
>> >>>>> >> >> > override default options (
>> >>>>> >> >> >
>> >>>>> >> >>
>> >>>>> >>
>> >>>>>
>> http://q4md-forcefieldtools.org/REDS-Development/Demo2-Files/Project.config
>> >>>>> >> >> ).
>> >>>>> >> >> > There is a portion of the Project.config file that deals with
>> >>>>> >> >> > intramolecular constraints (INTRA-MCC). My questions are:
>> >>>>> >> >>
>> >>>>> >> >> INTRA-MCC in R.E.D. Perl = INTRA-MCC1 in R.E.D. Python,
>> >>>>> which is the
>> >>>>> >> >> classical way to set up intra-molecular charge constraint...
>> >>>>> >> >> See Cieplak et al. in
>> >>>>> >> >>
>> >>>>>
>> http://www3.interscience.wiley.com/cgi-bin/abstract/109583237/ABSTRACT
>> >>>>> >> >>
>> >>>>> >> >> > - How does one decide what value to assign for
>> intramolecular
>> >>>>> >> >> > constraints?
>> >>>>> >> >>
>> >>>>> >> >> In general when one wants to remove a part of a whole molecule
>> to
>> >>>>> >> >> generate a molecular fragment the intra-mcc total charge
>> >>>>> value takes
>> >>>>> >> >> an integer value (& this integer value is very often zero).
>> >>>>> >> >>
>> >>>>> >> >> See also a discussion about the relation between the chemical
>> group
>> >>>>> >> >> involved in the constraint & the total charge of this
>> constraint:
>> >>>>> >> >> http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2918240/
>> >>>>> >> >>
>> >>>>> >> >> > - How does one decide whether or not to keep the
>> >>>>> intra-molecular
>> >>>>> >> >> charge
>> >>>>> >> >> > constraints through the second stage of fitting (i.e.,
>> >>>>> INTRA-MCC1
>> >>>>> >> vs.
>> >>>>> >> >> > INTRA-MCC2)?
>> >>>>> >> >>
>> >>>>> >> >> In general you want to use the way defined by Cieplak et al.
>> i.e.
>> >>>>> >> >> INTRA-MCC1
>> >>>>> >> >> See
>> >>>>> >>
>> http://www3.interscience.wiley.com/cgi-bin/abstract/109583237/ABSTRACT
>> >>>>> >> >>
>> >>>>> >> >> run both types of intra-mcc on the CH3CO group in a dipeptide
>> model
>> >>>>> to
>> >>>>> >> >> study the differences; in short an intra-mcc1 is only
>> >>>>> applied in the
>> >>>>> >> >> 1st resp input, while intra-mcc2 is applied in both resp
>> >>>>> inputs; the
>> >>>>> >> >> later allows charge equivalencing for CH2 and/or CH3 groups
>> when
>> >>>>> >> >> involved in the constraint...
>> >>>>> >> >>
>> >>>>> >> >> > If my questions are trivial ones, please direct me to the
>> >>>>> appropriate
>> >>>>> >> >> > references.
>> >>>>> >> >>
>> >>>>> >> >> no trivial question - do not worry ;-)
>> >>>>> >> >>
>> >>>>> >> >> R.E.D. Python is not yet published - still coding new
>> features...
>> >>>>> >> >> R.E.D. Perl is published:
>> >>>>> >> >> http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2918240/
>> >>>>> >> >>
>> >>>>> >> >> regards, Francois




_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Fri Dec 13 2013 - 08:00:02 PST
Custom Search