Re: [AMBER] CPPTRAJ, Total H-Bond Count Question

From: Daniel Roe via AMBER <amber.ambermd.org>
Date: Fri, 22 Jul 2022 15:41:37 -0400

OK - I've checked the first few frames and as far as I can tell
cpptraj and vmd are giving the same results for # hbonds calcd (I
modified your cpptraj script to give the number of hydrogen bonds vs
time via 'out numhbvtime.dat'):

$ head HBA100Kvmdhbonds.dat
0 386
1 381
2 370
3 362
4 345
5 371
6 365
7 372
8 355
9 362
$ awk '{print $1 "\t" $2;}' numhbvtime.dat
#Frame 100K[UU]
1 386
2 381
3 370
4 361
5 345
6 371
7 365
8 372
9 355
10 362

I also note you're writing out the average solvent-solute hbonds and
bridging info, but not the average solute-solute hbonds (easily
accomplished by adding 'avgout hbavg.dat').

As far as I can tell, both programs are producing equivalent results.
I'm guessing any differences you see is in how VMD groups the data
into the details (HBA100KvmdDetails.dat) file. The VMD code makes no
distinction between hydrogens, while cpptraj does. For example, if you
grab all the Lys 108 data out of the VMD file you get things like:

$ grep LYS108 HBA100KvmdDetails.dat
LYS108-Side-NZ ASP340-Side-OD2 33.28%
SER337-Main-N LYS108-Main-O 32.17%
LYS108-Side-NZ GLU104-Main-O 0.39%
LYS108-Side-NZ LYS91-Main-O 0.89%
LYS108-Side-NZ SER103-Main-O 0.17%

Whereas if you grab the Lys 108 data from cpptraj output you get:

$ grep LYS_108 hbavg.dat
LYS_108.O SER_337.H SER_337.N 579 0.3217
  2.8196 163.1371
ASP_340.OD2 LYS_108.HZ2 LYS_108.NZ 206 0.1144
  2.7326 166.1036
ASP_340.OD2 LYS_108.HZ3 LYS_108.NZ 205 0.1139
  2.7263 166.3880
ASP_340.OD2 LYS_108.HZ1 LYS_108.NZ 188 0.1044
  2.7293 165.7254
LYS_91.O LYS_108.HZ3 LYS_108.NZ 6 0.0033
  2.9093 165.3931
LYS_91.O LYS_108.HZ2 LYS_108.NZ 5 0.0028
  2.9101 163.7331
LYS_91.O LYS_108.HZ1 LYS_108.NZ 5 0.0028
  2.9107 167.9775
GLU_104.O LYS_108.HZ1 LYS_108.NZ 4 0.0022
  2.7859 161.7393
GLU_104.O LYS_108.HZ3 LYS_108.NZ 2 0.0011
  2.7666 160.9960
GLU_104.O LYS_108.HZ2 LYS_108.NZ 1 0.0006
  2.7385 166.7346
SER_103.O LYS_108.HZ3 LYS_108.NZ 1 0.0006
  2.7500 161.4854
SER_103.O LYS_108.HZ2 LYS_108.NZ 1 0.0006
  2.8366 162.1290
SER_103.O LYS_108.HZ1 LYS_108.NZ 1 0.0006
  2.9242 160.6470

Note that when there is only 1 hydrogen the results between the 2
codes are equivalent, but cpptraj breaks up the data per-hydrogen
while VMD does not. So it's just a difference in how the frames are
counted.

Hopefully this clears things up.

-Dan

On Fri, Jul 22, 2022 at 2:06 PM Daniel Roe <daniel.r.roe.gmail.com> wrote:
>
> Hi,
>
> It gave me an "Error downloading" before but seems to be working now,
> so the problem was likely a transient one on my end.
>
> -Dan
>
> On Fri, Jul 22, 2022 at 12:34 PM Matthew Guberman-Pfeffer
> <matthew.guberman-pfeffer.uconn.edu> wrote:
> >
> > Hi Dan, The link still works for me. Can you please let me know if it doesn’t work for you. What’s the message you see?
> > Best,
> > Matthew
> >
> > > On Jul 22, 2022, at 10:18 AM, Daniel Roe <daniel.r.roe.gmail.com> wrote:
> > >
> > > *Message sent from a system outside of UConn.*
> > >
> > >
> > > Hi,
> > >
> > > I don't think I got to the link in time. Any chance you can send a new
> > > one? Thanks,
> > >
> > > -Dan
> > >
> > > On Thu, Jul 7, 2022 at 4:13 PM Matthew Guberman-Pfeffer
> > > <matthew.guberman-pfeffer.uconn.edu> wrote:
> > >>
> > >> Hi Dan, I’m now looking at the details of the H-bonds instead of the total number, and I—to the best of my knowledge—see a difference between the output for H-bond analysis in Vmd and CPPTRAJ after setting the criteria to be equivalent. Below is a dropbox link to a directory with a ReadMe.txt file that explains what everything is and how I am comparing the outputs.
> > >>
> > >> Dropbox link: https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.dropbox.com%2Fsh%2Fghiwb6r4qqsl0k0%2FAAC-ZEARN3nQF-_HDtjKTU8ma%3Fdl%3D0&amp;data=05%7C01%7C%7C2ecf8ecf5c3846922db108da6bed04f5%7C17f1a87e2a254eaab9df9d439034b080%7C0%7C0%7C637940962970768008%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&amp;sdata=pxrdTUD3aM%2FzQ0c3UILgDjzg%2FiXUkRw4zpGLM0c329w%3D&amp;reserved=0
> > >>
> > >> Basically, I see that some H-bonds (all with an occupancy less than 10% but one as high as 7%) are found by VMD and not CPPTRAJ. I’m not saying either program is right/wrong; I only want to understand the origin of the difference.
> > >>
> > >> Best,
> > >> Matthew
> > >>
> > >>
> > >> On Jul 7, 2022, at 12:24 PM, Matthew Guberman-Pfeffer <matthew.guberman-pfeffer.uconn.edu> wrote:
> > >>
> > >> Dan, for you and others on the list, I want to revise my original post: If the angle and distance criteria are set to 160° and 3 Å, respectively, CPPTRAJ gives exactly the same H-bond count as VMD. I was mistaken before because I looked at the total number of H-bonds given by CPPTRAJ and the total number of unique H-bonds found by VMD. The first of these is on a per-frame basis, the latter is over the entire trajectory, which is why I thought VMD was finding many more H-bonds.
> > >>
> > >> Best,
> > >> Matthew
> > >>
> > >>
> > >> On Jun 24, 2022, at 3:36 PM, Daniel Roe <daniel.r.roe.gmail.com> wrote:
> > >>
> > >>
> > >> *Message sent from a system outside of UConn.*
> > >>
> > >>
> > >> Both cpptraj and vmd use simple geometric criteria - as long as the criteria match the results should be the same. Please send me off list a topology and coordinates so I can reproduce. Please also send the exact cpptraj input.
> > >>
> > >> On Fri, Jun 24, 2022 at 3:13 PM Matthew Guberman-Pfeffer <matthew.guberman-pfeffer.uconn.edu> wrote:
> > >>>
> > >>> Hi Dan,
> > >>>
> > >>> Yes, I noticed that CPPTRAJ and VMD use different definitions of the angles. I just worked in the reverse direction. I changed the angle in CPPTRAJ to be > 160° to match the VMD definition of < 20°. I think you were suggesting to change the VMD default to <45° to match the CPPTRAJ default of >135°.
> > >>>
> > >>> Using <20° in VMD and >160° in CPPTRAJ, I see a huge difference in the number of H-bonds (<300 by CPPTRAJ and almost 3000 by VMD). VMD also counts a huge number of H-bonds with water, whereas CPPTRAJ only reports a few hundred. Is something else besides the angle definition, which I think I already addressed, going on? Should I send you files so you can look at the comparison?
> > >>>
> > >>> Best,
> > >>> Matthew
> > >>>
> > >>>
> > >>>> On Jun 24, 2022, at 2:51 PM, Daniel Roe <daniel.r.roe.gmail.com> wrote:
> > >>>>
> > >>>> *Message sent from a system outside of UConn.*
> > >>>>
> > >>>>
> > >>>> Hi,
> > >>>>
> > >>>> On Mon, Jun 20, 2022 at 8:23 AM Matthew Guberman-Pfeffer via AMBER
> > >>>> <amber.ambermd.org> wrote:
> > >>>>>
> > >>>>> I’m trying to understand why the number is so small. If I run the H-bond analysis in VMD, I find ~2000 unique H-bonds. (I made sure to adjust the distance and angle criteria in CPPTRAJ and VMD to be the same for this comparison.)
> > >>>>
> > >>>> Make sure that the angle cutoff in particular is the same; cpptraj and
> > >>>> VMD defines the cutoffs differently. A valid hydrogen bond for cpptraj
> > >>>> has a D-H-A angle greater than the cutoff (default 135). In VMD it's
> > >>>> *deviation* from linear (180 deg.), so the cpptraj default of 135 is
> > >>>> 45 for VMD.
> > >>>>
> > >>>> Does that make sense?
> > >>>>
> > >>>> -Dan
> > >>>>
> > >>>>>
> > >>>>> Any insight would be much appreciated.
> > >>>>>
> > >>>>> Best,
> > >>>>> Matthew
> > >>>>>
> > >>>>>
> > >>>>> _______________________________________________
> > >>>>> AMBER mailing list
> > >>>>> AMBER.ambermd.org
> > >>>>> https://nam10.safelinks.protection.outlook.com/?url=http%3A%2F%2Flists.ambermd.org%2Fmailman%2Flistinfo%2Famber&amp;data=05%7C01%7C%7C2ecf8ecf5c3846922db108da6bed04f5%7C17f1a87e2a254eaab9df9d439034b080%7C0%7C0%7C637940962970768008%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&amp;sdata=hC7Ax1186qTCiJgR168aoAjzC84%2Bg32TBcpOrl7bObI%3D&amp;reserved=0
> > >>
> > >>
> > >>
> >

_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu Aug 04 2022 - 13:35:32 PDT
Custom Search