Hi Jason,
Thank you very much for the detailed answer. Since I am trying to use an
implicit solvent model, I computed the _ReferenceEnergy term using TI,
followed by typing in all the atom names and atomic charges into the
residues.py script accordingly:
# STP
refene1 = _ReferenceEnergy(igb2=0, igb5=0, igb8=0)
refene1.solvent_energies()
refene1.dielc2_energies()
refene1.dielc2.solvent_energies()
refene2 = _ReferenceEnergy(igb2=0, igb5=0, igb8=-83.6047)
refene2.solvent_energies()
refene2.dielc2_energies()
refene2.dielc2.solvent_energies()
refene2.set_pKa(5.2, deprotonated=True)
STP = TitratableResidue('STP', ['C2', 'H2', 'C1', 'H1', 'C5', 'H4', 'N1',
'H32', 'C4', 'H3', 'C3', 'C6', 'N3', 'C7', 'H5', 'C8',
'H6', 'N2', 'C9', 'N4', 'H7', 'C10', 'C15', 'H10',
'C11',
'C29', 'H29', 'H30', 'H31', 'C12', 'H8', 'C13', 'H9',
'C14', 'N5', 'H11', 'C16', 'O1', 'C17', 'C22', 'C21',
'H14', 'H15', 'C18', 'H12', 'C19', 'H13', 'C20',
'C23',
'H16', 'H17', 'N6', 'C24', 'H18', 'H19', 'C25', 'H20',
'H21', 'N7', 'C26', 'H22', 'H23', 'H24', 'C27', 'H25',
'H26', 'C28', 'H27', 'H28'], pka=5.2)
STP.add_state(protcnt=1, refene=refene1, # protonated
charges=[0.0280, 0.1765, -0.1875, 0.1944, 0.0576, 0.1992,
-0.2441,
0.3707, 0.1533, 0.1729, -0.0609, 0.6793, -0.9497, -0.8008,
0.2475,
0.5590, 0.0908, -0.7663, 1.1941, -0.7322, 0.4048, 0.1658,
-0.2871,
0.1961, 0.1087, -0.2446, 0.0813, 0.0813, 0.0813, -0.2257,
0.1857,
-0.3276, 0.2012, 0.3540, -0.6426, 0.3197, 0.8032, -0.5972,
-0.1946,
-0.1165, -0.1793, 0.1363, 0.1459, -0.1165, 0.1459, -0.1793,
0.1363,
0.1836, -0.2068, 0.1136, 0.1136, -0.1493, -0.2497, 0.1392,
0.1392,
-0.0799, 0.0993, 0.0993, -0.2388, -0.2242, 0.0984, 0.0984,
0.0984,
-0.0799, 0.0993, 0.0993, -0.2497, 0.1392, 0.1392])
STP.add_state(protcnt=0, refene=refene2, # deprotonated
charges=[0.1622, 0.1160, -0.5119, 0.1919, 0.4665, 0.0459,
-0.6915,
0.0000, 0.5046, 0.0454, -0.3953, 0.7772, -0.9263, -0.7760,
0.2223,
0.4890, 0.0827, -0.7788, 1.1570, -0.7288, 0.3966, 0.1834,
-0.2969,
0.2001, 0.1065, -0.2755, 0.0880, 0.0880, 0.0880, -0.2161,
0.1740,
-0.3502, 0.1925, 0.3551, -0.6481, 0.3136, 0.7933, -0.5879,
-0.1766,
-0.1168, -0.1824, 0.1327, 0.1471, -0.1168, 0.1471, -0.1824,
0.1327,
0.1640, -0.1339, 0.0881, 0.0881, -0.1668, -0.2420, 0.1386,
0.1386,
-0.1029, 0.1057, 0.1057, -0.2352, -0.2196, 0.0956, 0.0956,
0.0956,
-0.1029, 0.1057, 0.1057, -0.2420, 0.1386, 0.1386])
STP.check()
However, when I use cpinutil.py, the program says "CPIN generation
complete!", without any warning, while the generated cpin file doesn't have
any data, like this:
&CNSTPH
CHRGDAT=
PROTCNT=
RESNAME='System: Unknown',
RESSTATE=
STATENE=
TRESCNT=0,
/
I wonder what is wrong with my residues.py. Would you have any idea? Thank
you.
Best,
Olivia
On Wed, Aug 27, 2014 at 5:26 PM, Jason Swails <jason.swails.gmail.com>
wrote:
> On Wed, Aug 27, 2014 at 3:30 PM, M Olivia Kim <olivijuly.gmail.com> wrote:
>
> > Hello,
> >
> > I was trying to parameterize a nonstandard residue for constant pH MD in
> > Amber 14 and noticed that
> > $AMBERHOME/AmberTools/src/etc/cpinutils/residues.py script has been
> > primarily changed from the corresponding cpin_data.py in Amber 12 or
> older
> > version.
> >
> > Specifically, in the old cpin_data.py, I only needed to get the "relene"
> > term, which I obtain from TI computation. I guess this term corresponds
> to
> > the _ReferenceEnergy term in the residues.py script in Amber 14. However,
> > in the residues.py, there are a couple additional terms that I need to
> > define: solvent_energies and dielec2_energies.
> >
> > Would anyone be able to explain 1) what those terms mean; 2) how I can
> > obtain the values for each term; and 3) how they affect the total
> reference
> > energy? Thank you!
> >
>
> 1) solvent_energies is the reference energy for running CpHMD in explicit
> solvent. It is usually quite close to the "normal" TI in GB (since the
> protonation state moves are still done in GB). But since conformational
> ensembles may be slightly different, so too are the reference energies.
>
> dielec2_energies are the reference energies if you use an internal
> dielectric of 2. That was something I was trying when I couldn't get good
> results with CpHMD in explicit solvent. As the problem turned out to be
> the effective GB radii increased by the 4 dummy hydrogens on the
> carboxylate residues, I'd not recommend the dielectric 2 be used by
> anybody.
>
> 2) For the dielec2_energies, you can do a TI setting intdiel=2. If you
> aren't planning on using it though, there's no point. For the
> solvent_energies, you will need to use the original GB reference energies,
> run a titration (with pH-REMD for accelerated convergence), and then adjust
> the reference energy so that the computed pKa matches the reference pKa.
> This can be done in a single step (i.e., you don't need to iterate to get
> the correct reference energy until the pKa is properly predicted) -- I can
> describe exactly how if you need me to.
>
> In all cases, it's best to run a quick test titration (on GPUs you can do
> hundreds of nanoseconds on a model compound in a few minutes) and then
> adjust the reference energy to get the experimental pKa.
>
> 3) They are different reference energies -- if you tell cpinutil.py to give
> you a CPIN file that you plan to use in an explicit solvent environment, it
> uses the solvent_energies instead. Likewise if you plan to use an internal
> dielectric of 2. You only have to add the energies that you actually plan
> to use, though.
>
> Hope this helps,
> Jason
>
> --
> Jason M. Swails
> BioMaPS,
> Rutgers University
> Postdoctoral Researcher
> _______________________________________________
> 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 Sep 17 2014 - 11:30:03 PDT
