Hi Mark <div><br></div><div>It worked ! pdb2gmx is now properly processing the pdb, and moreover the potential energy from gmx is within 1.27 kJ/mol of that from Charmm. More on energy comparison shortly. </div><div><br></div>
<div>However, further tinkering of the rtp and pdb file was required to generate the topology: </div><div>1. </div><div>The original rtp entry which had HT1/2/3 for the methyl hydrogens was giving the error</div><div>"</div>
<div>Fatal error:</div><div>Atom H1 in residue ACE 1 was not found in rtp entry ACE with 6 atoms</div><div>while sorting atoms.</div><div> "</div><div>even though neither rtp nor pdb had H1 atom. I do not understand why it was insisting on H1/2/3. It worked upon changing HT1/2/3 to H1/2/3. </div>
<div><br></div><div>2. </div><div>I changed bond entry 'N CAT' -> 'N CT' in rtp entry of CT3. Guess CAT was a typo. </div><div><br></div><div>3. </div><div>The improper dihedral in CT3</div><div>' -C CT N -O ; commented to match with charmm psf '</div>
<div>was absent from the charmm psf. So I commented it. </div><div><br></div><div>fyi - final pdb and rtp entries are listed below. </div><div><br></div><div>Now, regarding the energies, I compared the different bonded and non-bonded energy components from charmm with those from gmx. </div>
<div>Bonds, dihedral and coulomb energies matched to within 0.01 kJ/mol which was great! </div><div>The discrepancy is in angles (U-B) where gmx value is higher by 1.6652 kJ/mol and LJ where gmx is lower by 0.459 kJ/mol. </div>
<div><br></div><div>Any idea on how to fix the angles (not too bothered by LJ) ? </div><div>I checked that the total number of angles in charmm psf and gmx top are the same (=54), but havn't yet compared individual entries. </div>
<div><br></div><div>btw, would you (or someone else here) have a sample charmm (inp, crd, psf) and gmx (mdp, gro, top) files where the two energies show a better match? It would help me with this troubleshooting. </div><div>
<br></div><div>Thanks for your help</div><div>Sandeep </div><div><br></div><div>================final rtp entries============================</div><div><div><div>[ ACE ]</div><div> [ atoms ]</div><div> CT3 CT3 -0.270 0</div>
<div> H1 HA 0.090 1</div><div> H2 HA 0.090 2</div><div> H3 HA 0.090 3</div><div> C C 0.510 4</div><div> O O -0.510 5</div><div> [ bonds ]</div>
<div> C CT3</div><div> C +N</div><div> CT3 H1</div><div> CT3 H2</div><div> CT3 H3</div><div> O C</div><div> [ impropers ]</div><div> C CT3 +N O</div>
<div><br></div></div></div><div><div>[ CT3 ]</div><div>; this can also be done with the .c.tdb, but the atom naming is different</div><div>; and this can matter</div><div> [ atoms ]</div><div> N NH1 -0.470 0</div>
<div> HN H 0.310 1</div><div> CT CT3 -0.110 2</div><div> H1 HA 0.090 3</div><div> H2 HA 0.090 4</div><div> H3 HA 0.090 5</div><div> [ bonds ]</div>
<div> -C N</div><div> N HN</div><div> N CT</div><div> CT H1</div><div> CT H2</div><div> CT H3</div><div> [ impropers ]</div><div> N -C CT HN</div>
<div>; -C CT N -O ; commented to match with charmm psf </div></div><div>==============final rtp entries==================</div><div><br></div><div>==============final pdb======================</div><div>
<div><div>REMARK NONE * </div><div>REMARK DATE: 8/10/11 5:30:23 CREATED BY USER: ss2029 </div><div>ATOM 1 CT3 ACE 1 -2.160 0.537 0.930 1.00 0.00 DIAL</div>
<div>ATOM 2 H1 ACE 1 -2.466 0.003 0.005 1.00 0.00 DIAL</div><div>ATOM 3 H2 ACE 1 -2.562 1.572 0.910 1.00 0.00 DIAL</div><div>ATOM 4 H3 ACE 1 -2.562 0.001 1.816 1.00 0.00 DIAL</div>
<div>ATOM 5 C ACE 1 -0.672 0.582 1.009 1.00 0.00 DIAL</div><div>ATOM 6 O ACE 1 -0.105 1.128 1.954 1.00 0.00 DIAL</div><div>ATOM 7 N ALA 2 0.000 0.000 0.000 1.00 0.00 DIAL</div>
<div>ATOM 8 HN ALA 2 -0.453 -0.444 -0.769 1.00 0.00 DIAL</div><div>ATOM 9 CA ALA 2 1.459 0.000 0.000 1.00 0.00 DIAL</div><div>ATOM 10 HA ALA 2 1.812 -0.497 0.897 1.00 0.00 DIAL</div>
<div>ATOM 11 CB ALA 2 1.949 -0.834 -1.207 1.00 0.00 DIAL</div><div>ATOM 12 HB1 ALA 2 1.514 -1.854 -1.155 1.00 0.00 DIAL</div><div>ATOM 13 HB2 ALA 2 1.628 -0.373 -2.167 1.00 0.00 DIAL</div>
<div>ATOM 14 HB3 ALA 2 3.056 -0.932 -1.214 1.00 0.00 DIAL</div><div>ATOM 15 C ALA 2 2.096 1.401 0.000 1.00 0.00 DIAL</div><div>ATOM 16 O ALA 2 1.425 2.432 0.000 1.00 0.00 DIAL</div>
<div>ATOM 17 N ALA 3 3.451 1.458 0.000 1.00 0.00 DIAL</div><div>ATOM 18 HN ALA 3 3.954 0.594 0.000 1.00 0.00 DIAL</div><div>ATOM 19 CA ALA 3 4.276 2.664 0.000 1.00 0.00 DIAL</div>
<div>ATOM 20 HA ALA 3 4.065 3.236 -0.897 1.00 0.00 DIAL</div><div>ATOM 21 CB ALA 3 3.864 3.539 1.207 1.00 0.00 DIAL</div><div>ATOM 22 HB1 ALA 3 2.776 3.756 1.155 1.00 0.00 DIAL</div>
<div>ATOM 23 HB2 ALA 3 4.063 3.014 2.167 1.00 0.00 DIAL</div><div>ATOM 24 HB3 ALA 3 4.407 4.508 1.214 1.00 0.00 DIAL</div><div>ATOM 25 C ALA 3 5.791 2.399 0.000 1.00 0.00 DIAL</div>
<div>ATOM 26 O ALA 3 6.597 3.328 0.000 1.00 0.00 DIAL</div><div>ATOM 27 N CT3 4 6.175 1.110 0.000 1.00 0.00 DIAL</div><div>ATOM 28 HN CT3 4 5.528 0.351 0.000 1.00 0.00 DIAL</div>
<div>ATOM 29 CT CT3 4 7.566 0.777 0.000 1.00 0.00 DIAL</div><div>ATOM 30 H1 CT3 4 8.048 1.201 0.907 1.00 0.00 DIAL</div><div><div>ATOM 31 H2 CT3 4 7.683 -0.327 0.000 1.00 0.00 DIAL</div>
<div>ATOM 32 H3 CT3 4 8.048 1.201 -0.907 1.00 0.00 DIAL</div><div>TER 33 CT3 4</div><div>END</div><div>==============final pdb======================</div></div></div></div><div><br>
</div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div> </div><div><br></div><div><br></div><div><br><br><div class="gmail_quote">On Fri, Sep 9, 2011 at 11:35 AM, Mark Abraham <span dir="ltr"><<a href="mailto:mark.abraham@anu.edu.au">mark.abraham@anu.edu.au</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;"><div class="im"><div> </div><div> </div><span>On 09/09/11, <b>Sandeep Somani </b><<a href="mailto:ssomani@gmail.com" target="_blank">ssomani@gmail.com</a>> wrote:</span>
<blockquote style="border-left:#00f 1px solid;padding-left:13px;margin-left:0px" type="cite">Hi Mark
<div><br></div><div>Thanks for the rtp definitions. I added the rtp defs to aminoacids.rtp and renamed residue and atom names in pdb file. </div><div><br></div><div>But pdb2gmx is still complaining about the terminal groups. It seems to have picked up ACE, but not CT3. </div>
<div>pdb2gmx error message and modified pdb file are listed below. Any idea ? Do I need to modify .c.tdb ? If so, do you have that entry as well :) ? </div></blockquote>
<div> </div></div><div>Ah, no. To get rid of the warning, you need to teach GROMACS that CT3 is a protein residue. (It knew about ACE from other force fields.) Copy residuetypes.dat from the share/top directory to your working directory, and add CT3 as Protein at the end.</div>
<div> </div><div>To get rid of the error, you need to tell pdb2gmx not to try to change your termini. Use pdb2gmx -ter, and choose "None" for both N and C termini.</div><div> </div><font color="#888888"><div>Mark </div>
</font><div><div></div><div class="h5"><blockquote style="border-left:#00f 1px solid;padding-left:13px;margin-left:0px" type="cite">
<div>
<div> </div><div><br></div><div>Thanks </div><div>Sandeep </div><div><br></div><div>=============pdb2gmx error=======================================</div><div>
<div>Back Off! I just backed up topol.top to ./#topol.top.3#</div><div>Processing chain 1 (32 atoms, 4 residues)</div><div>There are 3 donors and 3 acceptors</div><div>There are 4 hydrogen bonds</div><div>Identified residue ACE1 as a starting terminus.</div>
<div>Warning: Residue CT32 in chain has different type (Other) from starting residue ACE1 (Protein).</div><div>Identified residue ALA2 as a ending terminus.</div><div>8 out of 8 lines of specbond.dat converted successfully</div>
<div>Start terminus ACE-1: NH3+</div><div>End terminus ALA-2: COO-</div><div><br></div><div>-------------------------------------------------------</div><div>Program pdb2gmx_d, VERSION 4.5.4</div><div>Source code file: pdb2top.c, line: 1070</div>
<div><br></div><div>Fatal error:</div><div>atom N not found in buiding block 1ACE while combining tdb and rtp</div><div>For more information and tips for troubleshooting, please check the GROMACS</div><div>website at <a href="http://www.gromacs.org/Documentation/Errors" target="_blank">http://www.gromacs.org/Documentation/Errors</a></div>
<div>-------------------------------------------------------</div></div><div>================pdb2gmx error====================================</div><div><br></div><div>================modified pdb=====================================</div>
<div>
<div>
<div>REMARK NONE * </div><div>REMARK DATE: 8/10/11 5:30:23 CREATED BY USER: ss2029 </div><div><b>ATOM 1 CT3 ACE 1 -2.160 0.537 0.930 1.00 0.00 DIAL</b></div>
<div><b>ATOM 2 HT1 ACE 1 -2.466 0.003 0.005 1.00 0.00 DIAL</b></div><div><b>ATOM 3 HT2 ACE 1 -2.562 1.572 0.910 1.00 0.00 DIAL</b></div><div><b>ATOM 4 HT3 ACE 1 -2.562 0.001 1.816 1.00 0.00 DIAL</b></div>
<div><b>ATOM 5 C ACE 1 -0.672 0.582 1.009 1.00 0.00 DIAL</b></div><div><b>ATOM 6 O ACE 1 -0.105 1.128 1.954 1.00 0.00 DIAL</b></div><div>ATOM 7 N ALA 1 0.000 0.000 0.000 1.00 0.00 DIAL</div>
<div>ATOM 8 HN ALA 1 -0.453 -0.444 -0.769 1.00 0.00 DIAL</div><div>ATOM 9 CA ALA 1 1.459 0.000 0.000 1.00 0.00 DIAL</div><div>ATOM 10 HA ALA 1 1.812 -0.497 0.897 1.00 0.00 DIAL</div>
<div>ATOM 11 CB ALA 1 1.949 -0.834 -1.207 1.00 0.00 DIAL</div><div>ATOM 12 HB1 ALA 1 1.514 -1.854 -1.155 1.00 0.00 DIAL</div><div>ATOM 13 HB2 ALA 1 1.628 -0.373 -2.167 1.00 0.00 DIAL</div>
<div>ATOM 14 HB3 ALA 1 3.056 -0.932 -1.214 1.00 0.00 DIAL</div><div>ATOM 15 C ALA 1 2.096 1.401 0.000 1.00 0.00 DIAL</div><div>ATOM 16 O ALA 1 1.425 2.432 0.000 1.00 0.00 DIAL</div>
<div>ATOM 17 N ALA 2 3.451 1.458 0.000 1.00 0.00 DIAL</div><div>ATOM 18 HN ALA 2 3.954 0.594 0.000 1.00 0.00 DIAL</div><div>ATOM 19 CA ALA 2 4.276 2.664 0.000 1.00 0.00 DIAL</div>
<div>ATOM 20 HA ALA 2 4.065 3.236 -0.897 1.00 0.00 DIAL</div><div>ATOM 21 CB ALA 2 3.864 3.539 1.207 1.00 0.00 DIAL</div><div>ATOM 22 HB1 ALA 2 2.776 3.756 1.155 1.00 0.00 DIAL</div>
<div>ATOM 23 HB2 ALA 2 4.063 3.014 2.167 1.00 0.00 DIAL</div><div>ATOM 24 HB3 ALA 2 4.407 4.508 1.214 1.00 0.00 DIAL</div><div>ATOM 25 C ALA 2 5.791 2.399 0.000 1.00 0.00 DIAL</div>
<div>ATOM 26 O ALA 2 6.597 3.328 0.000 1.00 0.00 DIAL</div><div><b>ATOM 27 N CT3 2 6.175 1.110 0.000 1.00 0.00 DIAL</b></div><div><b>ATOM 28 HN CT3 2 5.528 0.351 0.000 1.00 0.00 DIAL</b></div>
<div><b>ATOM 29 CT CT3 2 7.566 0.777 0.000 1.00 0.00 DIAL</b></div><div><b>ATOM 30 HT1 CT3 2 8.048 1.201 0.907 1.00 0.00 DIAL</b></div></div><div>
<div><b>ATOM 31 HT2 CT3 2 7.683 -0.327 0.000 1.00 0.00 DIAL</b></div><div><b>ATOM 32 HT3 CT3 2 8.048 1.201 -0.907 1.00 0.00 DIAL</b></div><div>TER 33 CT3 2</div>
<div>END</div></div><div>================modified pdb=====================================</div></div><div><br></div><div><br></div><div><br><br>
<div class="gmail_quote">On Fri, Sep 9, 2011 at 2:52 AM, Mark Abraham <span dir="ltr"><<a href="mailto:Mark.Abraham@anu.edu.au" target="_blank">Mark.Abraham@anu.edu.au</a>></span> wrote:<br>
<blockquote style="border-left:#ccc 1px solid;margin:0px 0px 0px 0.8ex;padding-left:1ex" class="gmail_quote">
<div>On 9/09/2011 8:05 AM, Sandeep Somani wrote:<br>
<blockquote style="border-left:#ccc 1px solid;margin:0px 0px 0px 0.8ex;padding-left:1ex" class="gmail_quote">Hi<br><br>I am trying to set up a simulation for dialanine using charmm27 ff but am getting errors in pdb2gmx due to missing definitions in rtp file.<br>
<br>I created the pdb file (below) using Charmm for the sequence ACE-ALA-ALA-CT3.<br><br>pdb2gmx recognizes ALA atoms but not those of the terminal groups ACE and CT3.<br></blockquote><br></div>Yes, CHARMM in GROMACS has lacked these for some time. I asked over a year ago for them to be added, but that hasn't happened. You can add<br>
<br>[ ACE ]<br> [ atoms ]<br> CT3 CT3 -0.270 0<br> HT1 HA 0.090 1<br> HT2 HA 0.090 2<br> HT3 HA 0.090 3<br> C C 0.510 4<br> O O -0.510 5<br>
[ bonds ]<br> C CT3<br> C +N<br> CT3 HT31<br> CT3 HT32<br> CT3 HT33<br> O C<br> [ impropers ]<br> C CT3 +N O<br><br>[ CT3 ]<br>; this can also be done with the .c.tdb, but the atom naming is different<br>
; and this can matter<br> [ atoms ]<br> N NH1 -0.470 0<br> HN H 0.310 1<br> CT CT3 -0.110 2<br> HT1 HA 0.090 3<br> HT2 HA 0.090 4<br> HT3 HA 0.090 5<br>
[ bonds ]<br> -C N<br> N HN<br> N CAT<br> CT HT1<br> CT HT2<br> CT HT3<br><br> [ impropers ]<br> N -C CAT HN<br> -C CAT N -O<br>
<br>to the end of aminoacids.rtp in a local copy of the charmm27.ff folder to make this work.
<div><br><br>
<blockquote style="border-left:#ccc 1px solid;margin:0px 0px 0px 0.8ex;padding-left:1ex" class="gmail_quote"><br>Does someone have the appropriate definitions for these residues ?<br>( I have tried <a href="http://swissparam.ch/" target="_blank">http://swissparam.ch/</a> to generate an itp file for the molecule, but upon using it in GMX it does not give the same energy as that from charmm. )<br>
</blockquote><br></div>That will likely be a separate issue.<br><br>Mark<br><font color="#888888">-- <br>gmx-users mailing list <a href="mailto:gmx-users@gromacs.org" target="_blank">gmx-users@gromacs.org</a><br><a href="http://lists.gromacs.org/mailman/listinfo/gmx-users" target="_blank">http://lists.gromacs.org/<u></u>mailman/listinfo/gmx-users</a><br>
Please search the archive at <a href="http://www.gromacs.org/Support/Mailing_Lists/Search" target="_blank">http://www.gromacs.org/<u></u>Support/Mailing_Lists/Search</a> before posting!<br>Please don't post (un)subscribe requests to the list. Use the www interface or send it to <a href="mailto:gmx-users-request@gromacs.org" target="_blank">gmx-users-request@gromacs.org</a>.<br>
Can't post? Read <a href="http://www.gromacs.org/Support/Mailing_Lists" target="_blank">http://www.gromacs.org/<u></u>Support/Mailing_Lists</a><br></font></blockquote></div><br></div></div></blockquote>
<div> </div>
</div></div><br>--<br>
gmx-users mailing list <a href="mailto:gmx-users@gromacs.org">gmx-users@gromacs.org</a><br>
<a href="http://lists.gromacs.org/mailman/listinfo/gmx-users" target="_blank">http://lists.gromacs.org/mailman/listinfo/gmx-users</a><br>
Please search the archive at <a href="http://www.gromacs.org/Support/Mailing_Lists/Search" target="_blank">http://www.gromacs.org/Support/Mailing_Lists/Search</a> before posting!<br>
Please don't post (un)subscribe requests to the list. Use the<br>
www interface or send it to <a href="mailto:gmx-users-request@gromacs.org">gmx-users-request@gromacs.org</a>.<br>
Can't post? Read <a href="http://www.gromacs.org/Support/Mailing_Lists" target="_blank">http://www.gromacs.org/Support/Mailing_Lists</a><br></blockquote></div><br></div>