Re: [AMBER] Question for residues.py to parameterize a nonstandard residue for constant pH MD

From: Levi Pierce <levipierce.gmail.com>
Date: Sun, 28 Sep 2014 11:46:06 -0700

So it looks like this block in cpinutil.py is wonky??? We are hitting the
if statement and it gets evaluated like so

parm.ptr('natom') = 69
parm.parm_data['RESIDUE_POINTER'][resnum-1] = 1

so natoms = 68

Then natoms does not equal len(res.atom_list) and we never add it to the
list and we get an empty cpin file. If I comment out the continue it
generates the cpin file correctly what exactly is this if statement
checking for, should parm.parm_data['RESIDUE_POINTER'][resnum-1] actually
evaluate to 0 and not 1?

Thanks so much!

# Filter out termini (make sure the residue in the prmtop has as many
# atoms as the titratable residue defined in residues.py)

      if resnum == parm.ptr('nres'):
         natoms =
parm.ptr('natom')-parm.parm_data['RESIDUE_POINTER'][resnum-1]
         print "Checking", natoms, parm.ptr('natom'),
parm.parm_data['RESIDUE_POINTER'][resnum-1]
      else:

         natoms = (parm.parm_data['RESIDUE_POINTER'][resnum] -
                     parm.parm_data['RESIDUE_POINTER'][resnum-1])

      if natoms != len(res.atom_list):

          print "natoms final:", natoms, len(res.atom_list)
          continue
      # If we have gotten this far, add it to the list.
      main_reslist.add_residue(res, resnum,
                               parm.parm_data['RESIDUE_POINTER'][resnum-1])

      print "gets to adding"


On Wed, Sep 17, 2014 at 1:00 PM, M Olivia Kim <olivijuly.gmail.com> wrote:
>
> Hi,
>
> Yes I did add the new residue name to the titratable_residue list.
However,
> the problem still persists. Do you have any other suggestion or guess?
>
> Thanks,
> Olivia
>
> On Wed, Sep 17, 2014 at 2:40 PM, Jason Swails <jason.swails.gmail.com>
> wrote:
>
> > On Wed, 2014-09-17 at 14:11 -0400, M Olivia Kim wrote:
> > > 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.
> >
> > Don't forget to add the new residue name to the titratable_residues list
> > at the top of residues.py. That's the list that cpinutil reads from to
> > determine which residues are titratable.
> >
> > 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




--
--
Levi Pierce | Bioinformatics Engineer | Human Longevity, Inc. | 858-249-7559
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Sun Sep 28 2014 - 12:30:02 PDT
Custom Search