No template found for residue #2103

biosafetylvl5 · 2018-06-18T21:58:14Z

Hi! I'm trying to go a little bit farther than the getting started example and try out the same thing with a short oligonucleotide instead of a protein. However, I keep getting the same error of "No template found for residue."

I thought that this might be an issue with my PDB formation, so I tried a few oligos from the PDB (6GN4 and 6GPI) both of which have the same issue for me. I also created two new "oligos" - one which is single-stranded DNA and another which is just dTMP. Both give the same issue.

For the dTMP this is the error I get:

Traceback (most recent call last):
  File "C:\Users\alex\Documents\SSI\sims\sim1\sim1.py", line 9, in <module>
    nonbondedCutoff=1*nanometer, constraints=HBonds)
  File "C:\ProgramData\Anaconda3\lib\site-packages\simtk\openmm\app\forcefield.py", line 1128, in createSystem
    raise ValueError('No template found for residue %d (%s).  %s' % (res.index+1, res.name, _findMatchErrors(self, res)))
ValueError: No template found for residue 1 (DT).  This might mean your input topology is missing some atoms or bonds, or possibly that you are using the wrong force field.

which seems to be "recognising" the thymidine but still isn't able to go forward?

I'm sure I'm missing something small, but any help would be appreciated.

This is my code, almost exactly the same as the sample code:

from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
from sys import stdout

pdb = PDBFile('T.pdb')
forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')
system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME,
        nonbondedCutoff=1*nanometer, constraints=HBonds)
integrator = LangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
simulation = Simulation(pdb.topology, system, integrator)
simulation.context.setPositions(pdb.positions)
simulation.minimizeEnergy()
simulation.reporters.append(PDBReporter('output.pdb', 1000))
simulation.reporters.append(StateDataReporter(stdout, 1000, step=True, potentialEnergy=True, temperature=True))
simulation.step(10000)

The text was updated successfully, but these errors were encountered:

biosafetylvl5 · 2018-06-18T22:01:55Z

Here are the PDB file contents:
(these are the ones I made using Avogadro)

COMPND    UNNAMED
AUTHOR    GENERATED BY OPEN BABEL 2.3.90
ATOM      1 HTER DT  A   1      -1.547   9.377  -1.216  1.00  0.00           H
ATOM      2  OXT DT  A   1      -1.440   8.896  -0.401  1.00  0.00           O
ATOM      3  P   DT  A   1       0.000   8.910   0.000  1.00  0.00           P
ATOM      4  O1P DT  A   1       0.220  10.174   0.734  1.00  0.00           O1-
ATOM      5  O2P DT  A   1       0.790   8.732  -1.240  1.00  0.00           O
ATOM      6  O5' DT  A   1       0.250   7.669   0.971  1.00  0.00           O
ATOM      7  C5' DT  A   1      -0.690   7.422   2.046  1.00  0.00           C
ATOM      8  C4' DT  A   1       0.040   6.861   3.247  1.00  0.00           C
ATOM      9  O4' DT  A   1       0.250   5.428   3.036  1.00  0.00           O
ATOM     10  C3' DT  A   1       1.440   7.412   3.508  1.00  0.00           C
ATOM     11  O3' DT  A   1       1.830   7.271   4.868  1.00  0.00           O
ATOM     12  C2' DT  A   1       2.320   6.527   2.637  1.00  0.00           C
ATOM     13  C1' DT  A   1       1.610   5.184   2.732  1.00  0.00           C
ATOM     14  N1  DT  A   1       1.660   4.387   1.478  1.00  0.00           N
ATOM     15  C2  DT  A   1       1.810   3.022   1.600  1.00  0.00           C
ATOM     16  O2  DT  A   1       1.900   2.465   2.679  1.00  0.00           O
ATOM     17  N3  DT  A   1       1.850   2.324   0.409  1.00  0.00           N
ATOM     18  C4  DT  A   1       1.760   2.855  -0.855  1.00  0.00           C
ATOM     19  O4  DT  A   1       1.810   2.126  -1.853  1.00  0.00           O
ATOM     20  C5  DT  A   1       1.610   4.289  -0.887  1.00  0.00           C
ATOM     21  C5M DT  A   1       1.500   4.911  -2.246  1.00  0.00           C
ATOM     22  C6  DT  A   1       1.560   5.003   0.255  1.00  0.00           C
ATOM     23 H5'1 DT  A   1      -1.180   8.368   2.325  1.00  0.00           H
ATOM     24 H5'2 DT  A   1      -1.449   6.701   1.711  1.00  0.00           H
ATOM     25  H4' DT  A   1      -0.582   7.036   4.132  1.00  0.00           H
ATOM     26  H1' DT  A   1       2.073   4.598   3.539  1.00  0.00           H
ATOM     27 H2'1 DT  A   1       2.352   6.890   1.599  1.00  0.00           H
ATOM     28 H2'2 DT  A   1       3.344   6.467   3.037  1.00  0.00           H
ATOM     29  H3' DT  A   1       1.510   8.469   3.207  1.00  0.00           H
ATOM     30  H6  DT  A   1       1.439   6.085   0.207  1.00  0.00           H
ATOM     31  H51 DT  A   1       2.469   5.348  -2.531  1.00  0.00           H
ATOM     32  H52 DT  A   1       0.732   5.696  -2.229  1.00  0.00           H
ATOM     33  H53 DT  A   1       1.220   4.140  -2.980  1.00  0.00           H
ATOM     34  H3  DT  A   1       1.956   1.332   0.474  1.00  0.00           H
ATOM     35 HCAP DT  A   1       1.699   6.361   5.144  1.00  0.00           H  
CONECT    1    2
CONECT    2    1    3
CONECT    3    2    5    4    6
CONECT    4    3
CONECT    5    3
CONECT    6    3    7
CONECT    7    6   24   23    8
CONECT    8    7    9   10   25
CONECT    9   13    8
CONECT   10   12   29    8   11
CONECT   11   10   35
CONECT   12   27   13   28   10
CONECT   13   14   12    9   26
CONECT   14   22   15   13
CONECT   15   17   14   16
CONECT   16   15
CONECT   17   18   34   15
CONECT   18   19   20   17
CONECT   19   18
CONECT   20   21   18   22
CONECT   21   33   31   32   20
CONECT   22   20   30   14
CONECT   23    7
CONECT   24    7
CONECT   25    8
CONECT   26   13
CONECT   27   12
CONECT   28   12
CONECT   29   10
CONECT   30   22
CONECT   31   21
CONECT   32   21
CONECT   33   21
CONECT   34   17
CONECT   35   11
MASTER        0    0    0    0    0    0    0    0   35    0   35    0
END

And I've put the other one here as to not clutter this page

peastman · 2018-06-19T00:43:03Z

It looks like your nucleotide is terminated differently from what's defined in the force field. Amber14 requires nucleotides to be terminated in one specific way, and doesn't include parameters for any other form. Here is its definition of an isolated DT:

<Residue name="DTN">
  <Atom charge="0.4422" name="HO5'" type="DNA-HO"/>
  <Atom charge="-0.6318" name="O5'" type="DNA-OH"/>
  <Atom charge="-0.0069" name="C5'" type="DNA-CJ"/>
  <Atom charge="0.0754" name="H5'" type="DNA-H1"/>
  <Atom charge="0.0754" name="H5''" type="DNA-H1"/>
  <Atom charge="0.1629" name="C4'" type="DNA-CT"/>
  <Atom charge="0.1176" name="H4'" type="DNA-H1"/>
  <Atom charge="-0.3691" name="O4'" type="DNA-OS"/>
  <Atom charge="0.068" name="C1'" type="DNA-CT"/>
  <Atom charge="0.1804" name="H1'" type="DNA-H2"/>
  <Atom charge="-0.0239" name="N1" type="DNA-N*"/>
  <Atom charge="-0.2209" name="C6" type="DNA-CM"/>
  <Atom charge="0.2607" name="H6" type="DNA-H4"/>
  <Atom charge="0.0025" name="C5" type="DNA-CM"/>
  <Atom charge="-0.2269" name="C7" type="DNA-CT"/>
  <Atom charge="0.077" name="H71" type="DNA-HC"/>
  <Atom charge="0.077" name="H72" type="DNA-HC"/>
  <Atom charge="0.077" name="H73" type="DNA-HC"/>
  <Atom charge="0.5194" name="C4" type="DNA-C"/>
  <Atom charge="-0.5563" name="O4" type="DNA-O"/>
  <Atom charge="-0.434" name="N3" type="DNA-NA"/>
  <Atom charge="0.342" name="H3" type="DNA-H"/>
  <Atom charge="0.5677" name="C2" type="DNA-C"/>
  <Atom charge="-0.5881" name="O2" type="DNA-O"/>
  <Atom charge="0.0713" name="C3'" type="DNA-C7"/>
  <Atom charge="0.0985" name="H3'" type="DNA-H1"/>
  <Atom charge="-0.0854" name="C2'" type="DNA-CT"/>
  <Atom charge="0.0718" name="H2'" type="DNA-HC"/>
  <Atom charge="0.0718" name="H2''" type="DNA-HC"/>
  <Atom charge="-0.6549" name="O3'" type="DNA-OH"/>
  <Atom charge="0.4396" name="HO3'" type="DNA-HO"/>
  <Bond atomName1="HO5'" atomName2="O5'"/>
  <Bond atomName1="O5'" atomName2="C5'"/>
  <Bond atomName1="C5'" atomName2="H5'"/>
  <Bond atomName1="C5'" atomName2="H5''"/>
  <Bond atomName1="C5'" atomName2="C4'"/>
  <Bond atomName1="C4'" atomName2="H4'"/>
  <Bond atomName1="C4'" atomName2="O4'"/>
  <Bond atomName1="C4'" atomName2="C3'"/>
  <Bond atomName1="O4'" atomName2="C1'"/>
  <Bond atomName1="C1'" atomName2="H1'"/>
  <Bond atomName1="C1'" atomName2="N1"/>
  <Bond atomName1="C1'" atomName2="C2'"/>
  <Bond atomName1="N1" atomName2="C6"/>
  <Bond atomName1="N1" atomName2="C2"/>
  <Bond atomName1="C6" atomName2="H6"/>
  <Bond atomName1="C6" atomName2="C5"/>
  <Bond atomName1="C5" atomName2="C7"/>
  <Bond atomName1="C5" atomName2="C4"/>
  <Bond atomName1="C7" atomName2="H71"/>
  <Bond atomName1="C7" atomName2="H72"/>
  <Bond atomName1="C7" atomName2="H73"/>
  <Bond atomName1="C4" atomName2="O4"/>
  <Bond atomName1="C4" atomName2="N3"/>
  <Bond atomName1="N3" atomName2="H3"/>
  <Bond atomName1="N3" atomName2="C2"/>
  <Bond atomName1="C2" atomName2="O2"/>
  <Bond atomName1="C3'" atomName2="H3'"/>
  <Bond atomName1="C3'" atomName2="C2'"/>
  <Bond atomName1="C3'" atomName2="O3'"/>
  <Bond atomName1="C2'" atomName2="H2'"/>
  <Bond atomName1="C2'" atomName2="H2''"/>
  <Bond atomName1="O3'" atomName2="HO3'"/>
</Residue>

biosafetylvl5 · 2018-06-19T21:55:14Z

Ah, I see it now - having drawn things out - I had a 5' phosphate but Amber wanted a 5' hydroxyl group. I've included a PDB for an isolated DT and a short oligo for others who want an example/something easy to start out with. Thanks for the help!

PDBs in zip file can be found here

asgharrazavi · 2020-09-21T19:37:23Z

I am also struggling with terminal residue patching when using charmm forcefiled. Could you please comment on the difference between e.g. GLY, GLYN, and GLYC in the charmm36.xml file? Are the "GLYN" and "GLYC" terminal residue capping for GLY? How does the AllowPatch work?

<Residue name="GLY">
  <Atom charge="-0.47" name="N" type="NH1"/>
  <Atom charge="0.31" name="HN" type="H"/>
  <Atom charge="-0.02" name="CA" type="CT2"/>
  <Atom charge="0.09" name="HA1" type="HB2"/>
  <Atom charge="0.09" name="HA2" type="HB2"/>
  <Atom charge="0.51" name="C" type="C"/>
  <Atom charge="-0.51" name="O" type="O"/>
  <Bond atomName1="N" atomName2="HN"/>
  <Bond atomName1="N" atomName2="CA"/>
  <Bond atomName1="C" atomName2="CA"/>
  <Bond atomName1="CA" atomName2="HA1"/>
  <Bond atomName1="CA" atomName2="HA2"/>
  <Bond atomName1="O" atomName2="C"/>
  <ExternalBond atomName="N"/>
  <ExternalBond atomName="C"/>
  <AllowPatch name="GLYP"/>
  <AllowPatch name="ACE"/>
  <AllowPatch name="NGNE"/>
  <AllowPatch name="CTER"/>
  <AllowPatch name="CNEU"/>
  <AllowPatch name="CT2"/>
  <AllowPatch name="CT3"/>
</Residue>
<Residue name="GLYN">
  <Atom charge="-0.96" name="N" type="NG321"/>
  <Atom charge="0.34" name="HT1" type="HGPAM2"/>
  <Atom charge="0.34" name="HT2" type="HGPAM2"/>
  <Atom charge="0.18" name="CA" type="CG321"/>
  <Atom charge="0.05" name="HA1" type="HGA2"/>
  <Atom charge="0.05" name="HA2" type="HGA2"/>
  <Atom charge="0.75" name="C" type="CG2O2"/>
  <Atom charge="-0.55" name="OT1" type="OG2D1"/>
  <Atom charge="-0.63" name="OT2" type="OG311"/>
  <Atom charge="0.43" name="HO2" type="HGP1"/>
  <Bond atomName1="N" atomName2="CA"/>
  <Bond atomName1="CA" atomName2="C"/>
  <Bond atomName1="C" atomName2="OT2"/>
  <Bond atomName1="OT2" atomName2="HO2"/>
  <Bond atomName1="N" atomName2="HT1"/>
  <Bond atomName1="N" atomName2="HT2"/>
  <Bond atomName1="CA" atomName2="HA1"/>
  <Bond atomName1="CA" atomName2="HA2"/>
  <Bond atomName1="OT1" atomName2="C"/>
  <AllowPatch name="ACE"/>
</Residue>
<Residue name="GLYC">
  <Atom charge="-0.27" name="C3" type="CG331"/>
  <Atom charge="0.09" name="H31" type="HGA3"/>
  <Atom charge="0.09" name="H32" type="HGA3"/>
  <Atom charge="0.09" name="H33" type="HGA3"/>
  <Atom charge="0.08" name="C1" type="CG321"/>
  <Atom charge="0.09" name="H11" type="HGA2"/>
  <Atom charge="0.09" name="H12" type="HGA2"/>
  <Atom charge="-0.49" name="O4" type="OG302"/>
  <Atom charge="0.9" name="C5" type="CG2O2"/>
  <Atom charge="-0.63" name="O6" type="OG2D1"/>
  <Atom charge="-0.31" name="C7" type="CG331"/>
  <Atom charge="0.09" name="H71" type="HGA3"/>
  <Atom charge="0.09" name="H72" type="HGA3"/>
  <Atom charge="0.09" name="H73" type="HGA3"/>
  <Atom charge="0.17" name="C2" type="CG311"/>
  <Atom charge="0.09" name="H21" type="HGA1"/>
  <Atom charge="-0.49" name="O8" type="OG302"/>
  <Atom charge="0.9" name="C9" type="CG2O2"/>
  <Atom charge="-0.63" name="O10" type="OG2D1"/>
  <Atom charge="-0.31" name="C11" type="CG331"/>
  <Atom charge="0.09" name="H111" type="HGA3"/>
  <Atom charge="0.09" name="H112" type="HGA3"/>
  <Atom charge="0.09" name="H113" type="HGA3"/>
  <Bond atomName1="C1" atomName2="C2"/>
  <Bond atomName1="C2" atomName2="C3"/>
  <Bond atomName1="C1" atomName2="O4"/>
  <Bond atomName1="O4" atomName2="C5"/>
  <Bond atomName1="C5" atomName2="C7"/>
  <Bond atomName1="C5" atomName2="O6"/>
  <Bond atomName1="C2" atomName2="O8"/>
  <Bond atomName1="O8" atomName2="C9"/>
  <Bond atomName1="C9" atomName2="C11"/>
  <Bond atomName1="C9" atomName2="O10"/>
  <Bond atomName1="C1" atomName2="H11"/>
  <Bond atomName1="C1" atomName2="H12"/>
  <Bond atomName1="C2" atomName2="H21"/>
  <Bond atomName1="C3" atomName2="H31"/>
  <Bond atomName1="C3" atomName2="H32"/>
  <Bond atomName1="C3" atomName2="H33"/>
  <Bond atomName1="C7" atomName2="H71"/>
  <Bond atomName1="C7" atomName2="H72"/>
  <Bond atomName1="C7" atomName2="H73"/>
  <Bond atomName1="C11" atomName2="H111"/>
  <Bond atomName1="C11" atomName2="H112"/>
  <Bond atomName1="C11" atomName2="H113"/>
  <AllowPatch name="MET1"/>
  <AllowPatch name="MET2"/>
</Residue>

I have CONNECT information in my pdb file but still get errors about not being able to match residue to template. Here is an example error:

ValueError: No template found for residue 1 (MET).  The set of atoms is similar to AOBT, but it is missing 2 atoms.

Here is my MET residue:

ATOM      1  N   MET A   1      17.541  -6.380 -13.937  1.00 20.33           N
ANISOU    1  N   MET A   1     1924   3195   2604    303    186   -292
ATOM      2  CA  MET A   1      16.391  -5.472 -13.960  1.00 18.20           C
ANISOU    2  CA  MET A   1     1702   2890   2322    245    190   -254
ATOM      3  C   MET A   1      15.091  -6.280 -13.875  1.00 15.92           C
ANISOU    3  C   MET A   1     1493   2543   2013    263    152   -267
ATOM      4  O   MET A   1      14.875  -7.180 -14.676  1.00 14.01           O
ANISOU    4  O   MET A   1     1276   2306   1741    309    153   -296
ATOM      5  CB  MET A   1      16.426  -4.612 -15.223  1.00 21.42           C
ANISOU    5  CB  MET A   1     2097   3346   2697    223    253   -227
ATOM      6  CG  MET A   1      15.486  -3.433 -15.193  1.00 25.49           C
ANISOU    6  CG  MET A   1     2649   3831   3205    162    259   -183
ATOM      7  SD  MET A   1      16.161  -1.963 -16.023  1.00 31.29           S
ANISOU    7  SD  MET A   1     3339   4618   3934    111    331   -136
ATOM      8  CE  MET A   1      16.121  -2.465 -17.605  1.00 28.17           C
ANISOU    8  CE  MET A   1     2958   4271   3473    159    378   -147
ATOM      9  C1  MET A   1      18.820  -5.876 -14.053  1.00  0.00           C
ATOM     10  O1  MET A   1      19.774  -6.618 -14.035  1.00  0.00           O
ATOM     11  C2  MET A   1      19.017  -4.356 -14.205  1.00  0.00           C
ATOM     12  H   MET A   1      17.397  -7.375 -13.837  1.00  0.00           H
ATOM     13  HA  MET A   1      16.449  -4.816 -13.092  1.00  0.00           H
ATOM     14  HB3 MET A   1      16.199  -5.233 -16.090  1.00  0.00           H
ATOM     15  HB2 MET A   1      17.443  -4.259 -15.392  1.00  0.00           H
ATOM     16  HG3 MET A   1      15.247  -3.187 -14.159  1.00  0.00           H
ATOM     17  HG2 MET A   1      14.541  -3.711 -15.661  1.00  0.00           H
ATOM     18  HE1 MET A   1      16.505  -1.673 -18.248  1.00  0.00           H
ATOM     19  HE2 MET A   1      16.738  -3.356 -17.720  1.00  0.00           H
ATOM     20  HE3 MET A   1      15.094  -2.696 -17.888  1.00  0.00           H
ATOM     21  H1  MET A   1      18.046  -3.861 -14.199  1.00  0.00           H
ATOM     22  H2  MET A   1      19.620  -3.983 -13.378  1.00  0.00           H
ATOM     23  H3  MET A   1      19.525  -4.148 -15.147  1.00  0.00           H

peastman · 2020-09-22T17:09:33Z

GLY, GLYC, and GLYN are the versions for middle of a chain, C-terminal, and N-terminal respectively. GLY is unusual in this respect. Most amino acids just have one version and use patches for the terminal residues.

Your MET residue looks strange. What are the C1, C2, and O1 atoms? Those aren't part of a standard MET residue.

asgharrazavi · 2020-09-22T18:39:38Z

Sorry for reopening this issue, but it turns out I still have the patching issue. To make it more clear and simple, how can I add terminal patching to this simple 3 residue peptide?

CRYST1   41.604  116.842   91.486  90.00  90.00  90.00 C 2 2 21      8
MODEL        1
ATOM      1  CH3 ACE A   0      35.416   9.009 -18.663  1.00  0.00           C  
ATOM      2  C   ACE A   0      34.218   9.958 -18.693  1.00  0.00           C  
ATOM      3  O   ACE A   0      33.432   9.978 -17.748  1.00  0.00           O  
ATOM      4 1H   ACE A   0      35.403   8.426 -17.715  1.00  0.00           H  
ATOM      5 2H   ACE A   0      35.362   8.313 -19.530  1.00  0.00           H  
ATOM      6 3H   ACE A   0      36.359   9.596 -18.727  1.00  0.00           H  
ATOM      7  N   MET A   1      34.018  10.791 -19.775  1.00 20.33           N  
ATOM      8  CA  MET A   1      32.868  11.699 -19.798  1.00 18.20           C  
ATOM      9  C   MET A   1      31.568  10.891 -19.713  1.00 15.92           C  
ATOM     10  O   MET A   1      31.352   9.991 -20.514  1.00 14.01           O  
ATOM     11  CB  MET A   1      32.903  12.559 -21.061  1.00 21.42           C  
ATOM     12  CG  MET A   1      31.963  13.738 -21.031  1.00 25.49           C  
ATOM     13  SD  MET A   1      32.638  15.208 -21.861  1.00 31.29           S  
ATOM     14  CE  MET A   1      32.598  14.706 -23.443  1.00 28.17           C  
ATOM     15  H1  MET A   1      34.664  10.775 -20.551  1.00  0.00           H  
ATOM     16  HA  MET A   1      32.926  12.355 -18.930  1.00  0.00           H  
ATOM     17  HB3 MET A   1      32.676  11.938 -21.927  1.00  0.00           H  
ATOM     18  HB2 MET A   1      33.920  12.912 -21.229  1.00  0.00           H  
ATOM     19  HG3 MET A   1      31.724  13.984 -19.996  1.00  0.00           H  
ATOM     20  HG2 MET A   1      31.018  13.460 -21.499  1.00  0.00           H  
ATOM     21  HE1 MET A   1      32.982  15.498 -24.086  1.00  0.00           H  
ATOM     22  HE2 MET A   1      33.215  13.815 -23.558  1.00  0.00           H  
ATOM     23  HE3 MET A   1      31.571  14.475 -23.726  1.00  0.00           H  
ATOM     24  N   THR A   2      30.712  11.223 -18.741  1.00  9.98           N  
ATOM     25  CA  THR A   2      29.469  10.501 -18.503  1.00  8.85           C  
ATOM     26  C   THR A   2      28.342  11.463 -18.217  1.00 11.00           C  
ATOM     27  O   THR A   2      28.506  12.405 -17.424  1.00  9.24           O  
ATOM     28  CB  THR A   2      29.651   9.536 -17.306  1.00 15.81           C  
ATOM     29  OG1 THR A   2      30.789   8.703 -17.531  1.00 13.83           O  
ATOM     30  CG2 THR A   2      28.415   8.664 -17.034  1.00 15.49           C  
ATOM     31  H   THR A   2      30.915  12.004 -18.134  1.00  0.00           H  
ATOM     32  HA  THR A   2      29.222   9.920 -19.391  1.00  0.00           H  
ATOM     33  HB  THR A   2      29.843  10.135 -16.416  1.00  0.00           H  
ATOM     34  HG1 THR A   2      30.901   8.105 -16.788  1.00  0.00           H  
ATOM     35 HG21 THR A   2      28.610   8.012 -16.183  1.00  0.00           H  
ATOM     36 HG22 THR A   2      27.560   9.303 -16.813  1.00  0.00           H  
ATOM     37 HG23 THR A   2      28.197   8.058 -17.913  1.00  0.00           H  
ATOM     38  N   GLU A   3      27.186  11.203 -18.842  1.00  8.87           N  
ATOM     39  CA  GLU A   3      25.982  11.973 -18.580  1.00  7.86           C  
ATOM     40  C   GLU A   3      25.238  11.354 -17.409  1.00 11.09           C  
ATOM     41  O   GLU A   3      25.183  10.124 -17.284  1.00 10.54           O  
ATOM     42  CB  GLU A   3      25.055  11.978 -19.791  1.00  9.04           C  
ATOM     43  CG  GLU A   3      25.462  12.978 -20.839  1.00 12.16           C  
ATOM     44  CD  GLU A   3      24.509  13.090 -22.010  1.00 19.31           C  
ATOM     45  OE1 GLU A   3      24.886  13.774 -22.981  1.00 11.19           O  
ATOM     46  OE2 GLU A   3      23.403  12.500 -21.976  1.00 12.96           O1-
ATOM     47  H   GLU A   3      27.134  10.454 -19.517  1.00  0.00           H  
ATOM     48  HA  GLU A   3      26.256  12.998 -18.332  1.00  0.00           H  
ATOM     49  HB3 GLU A   3      24.036  12.188 -19.466  1.00  0.00           H  
ATOM     50  HB2 GLU A   3      25.030  10.981 -20.232  1.00  0.00           H  
ATOM     51  HG3 GLU A   3      26.458  12.730 -21.207  1.00  0.00           H  
ATOM     52  HG2 GLU A   3      25.581  13.958 -20.376  1.00  0.00           H  
ATOM     53  N   NMA A   3A     24.622  12.170 -16.482  1.00  0.00           N  
ATOM     54  CA  NMA A   3A     23.907  11.584 -15.358  1.00  0.00           C  
ATOM     55  H   NMA A   3A     24.667  13.173 -16.584  1.00  0.00           H  
ATOM     56 1HA  NMA A   3A     23.976  10.475 -15.428  1.00  0.00           H  
ATOM     57 2HA  NMA A   3A     24.361  11.929 -14.402  1.00  0.00           H  
ATOM     58 3HA  NMA A   3A     22.839  11.897 -15.386  1.00  0.00           H  
CONECT    1    4    5    6
CONECT    4    1
CONECT    5    1
CONECT    6    1
ENDMDL
END

If I use the above PDB file, I get this error:

raise ValueError('No template found for residue %d (%s).  %s' % (res.index+1, res.name, _findMatchErrors(self, res)))
ValueError: No template found for residue 1 (ACE).  The set of atoms is similar to ACET, but it is missing 1 atoms.

Just to be clear, I get the above error only when using chamm36.xml. It will work fine when using amber99sbildn.xml.

peastman · 2020-09-22T18:56:32Z

CHARMM handles terminal groups differently from Amber. Amber considers ACE to be a separate residue, while CHARMM treats it as a patch on the existing residue. So change those lines of the PDB file to say it's part of the MET residue.

asgharrazavi · 2020-09-22T19:07:39Z

Unfortunately, that also doesn't work:

PDB:

CRYST1   41.604  116.842   91.486  90.00  90.00  90.00 C 2 2 21      8
MODEL        1
ATOM      1  CH3 MET A   1      35.416   9.009 -18.663  1.00  0.00           C  
ATOM      2  C   MET A   1      34.218   9.958 -18.693  1.00  0.00           C  
ATOM      3  O   MET A   1      33.432   9.978 -17.748  1.00  0.00           O  
ATOM      4 1H   MET A   1      35.403   8.426 -17.715  1.00  0.00           H  
ATOM      5 2H   MET A   1      35.362   8.313 -19.530  1.00  0.00           H  
ATOM      6 3H   MET A   1      36.359   9.596 -18.727  1.00  0.00           H  
ATOM      7  N   MET A   1      34.018  10.791 -19.775  1.00 20.33           N  
ATOM      8  CA  MET A   1      32.868  11.699 -19.798  1.00 18.20           C  
ATOM      9  C   MET A   1      31.568  10.891 -19.713  1.00 15.92           C  
ATOM     10  O   MET A   1      31.352   9.991 -20.514  1.00 14.01           O  
ATOM     11  CB  MET A   1      32.903  12.559 -21.061  1.00 21.42           C  
ATOM     12  CG  MET A   1      31.963  13.738 -21.031  1.00 25.49           C  
ATOM     13  SD  MET A   1      32.638  15.208 -21.861  1.00 31.29           S  
ATOM     14  CE  MET A   1      32.598  14.706 -23.443  1.00 28.17           C  
ATOM     15  H1  MET A   1      34.664  10.775 -20.551  1.00  0.00           H  
ATOM     16  HA  MET A   1      32.926  12.355 -18.930  1.00  0.00           H  
ATOM     17  HB3 MET A   1      32.676  11.938 -21.927  1.00  0.00           H  
ATOM     18  HB2 MET A   1      33.920  12.912 -21.229  1.00  0.00           H  
ATOM     19  HG3 MET A   1      31.724  13.984 -19.996  1.00  0.00           H  
ATOM     20  HG2 MET A   1      31.018  13.460 -21.499  1.00  0.00           H  
ATOM     21  HE1 MET A   1      32.982  15.498 -24.086  1.00  0.00           H  
ATOM     22  HE2 MET A   1      33.215  13.815 -23.558  1.00  0.00           H  
ATOM     23  HE3 MET A   1      31.571  14.475 -23.726  1.00  0.00           H  
ATOM     24  N   THR A   2      30.712  11.223 -18.741  1.00  9.98           N  
ATOM     25  CA  THR A   2      29.469  10.501 -18.503  1.00  8.85           C  
ATOM     26  C   THR A   2      28.342  11.463 -18.217  1.00 11.00           C  
ATOM     27  O   THR A   2      28.506  12.405 -17.424  1.00  9.24           O  
ATOM     28  CB  THR A   2      29.651   9.536 -17.306  1.00 15.81           C  
ATOM     29  OG1 THR A   2      30.789   8.703 -17.531  1.00 13.83           O  
ATOM     30  CG2 THR A   2      28.415   8.664 -17.034  1.00 15.49           C  
ATOM     31  H   THR A   2      30.915  12.004 -18.134  1.00  0.00           H  
ATOM     32  HA  THR A   2      29.222   9.920 -19.391  1.00  0.00           H  
ATOM     33  HB  THR A   2      29.843  10.135 -16.416  1.00  0.00           H  
ATOM     34  HG1 THR A   2      30.901   8.105 -16.788  1.00  0.00           H  
ATOM     35 HG21 THR A   2      28.610   8.012 -16.183  1.00  0.00           H  
ATOM     36 HG22 THR A   2      27.560   9.303 -16.813  1.00  0.00           H  
ATOM     37 HG23 THR A   2      28.197   8.058 -17.913  1.00  0.00           H  
ATOM     38  N   GLU A   3      27.186  11.203 -18.842  1.00  8.87           N  
ATOM     39  CA  GLU A   3      25.982  11.973 -18.580  1.00  7.86           C  
ATOM     40  C   GLU A   3      25.238  11.354 -17.409  1.00 11.09           C  
ATOM     41  O   GLU A   3      25.183  10.124 -17.284  1.00 10.54           O  
ATOM     42  CB  GLU A   3      25.055  11.978 -19.791  1.00  9.04           C  
ATOM     43  CG  GLU A   3      25.462  12.978 -20.839  1.00 12.16           C  
ATOM     44  CD  GLU A   3      24.509  13.090 -22.010  1.00 19.31           C  
ATOM     45  OE1 GLU A   3      24.886  13.774 -22.981  1.00 11.19           O  
ATOM     46  OE2 GLU A   3      23.403  12.500 -21.976  1.00 12.96           O1-
ATOM     47  H   GLU A   3      27.134  10.454 -19.517  1.00  0.00           H  
ATOM     48  HA  GLU A   3      26.256  12.998 -18.332  1.00  0.00           H  
ATOM     49  HB3 GLU A   3      24.036  12.188 -19.466  1.00  0.00           H  
ATOM     50  HB2 GLU A   3      25.030  10.981 -20.232  1.00  0.00           H  
ATOM     51  HG3 GLU A   3      26.458  12.730 -21.207  1.00  0.00           H  
ATOM     52  HG2 GLU A   3      25.581  13.958 -20.376  1.00  0.00           H  
ATOM     53  N   NMA A   3A     24.622  12.170 -16.482  1.00  0.00           N  
ATOM     54  CA  NMA A   3A     23.907  11.584 -15.358  1.00  0.00           C  
ATOM     55  H   NMA A   3A     24.667  13.173 -16.584  1.00  0.00           H  
ATOM     56 1HA  NMA A   3A     23.976  10.475 -15.428  1.00  0.00           H  
ATOM     57 2HA  NMA A   3A     24.361  11.929 -14.402  1.00  0.00           H  
ATOM     58 3HA  NMA A   3A     22.839  11.897 -15.386  1.00  0.00           H  
CONECT    1    4    5    6
CONECT    4    1
CONECT    5    1
CONECT    6    1
ENDMDL
END

the code:

from simtk.openmm.app import PDBFile
from simtk.openmm import app

print('building protein and solvating...')

# Load the AMBER protein force field parameters through OpenMM.
omm_forcefield = app.ForceField('charmm36.xml', 'charmm36/water.xml')

# Load the solvated receptor from a PDB file.
prot_pdb_file = PDBFile('3res.pdb') 
modeller = app.Modeller(prot_pdb_file.topology, prot_pdb_file.positions)

print('Adding hydrogens...')
modeller.addHydrogens(omm_forcefield)

And the entire output:

building protein and solvating...
/home/asgharrazavi/miniconda3/lib/python3.7/site-packages/simtk/openmm/app/internal/pdbstructure.py:537: UserWarning: WARNING: duplicate atom (ATOM      9  C   MET A   1      31.568  10.891 -19.713  1.00 15.92           C  , ATOM      2  C   MET A   1      34.218   9.958 -18.693  1.00  0.00           C  )
  warnings.warn("WARNING: duplicate atom (%s, %s)" % (atom, old_atom._pdb_string(old_atom.serial_number, atom.alternate_location_indicator)))
/home/asgharrazavi/miniconda3/lib/python3.7/site-packages/simtk/openmm/app/internal/pdbstructure.py:537: UserWarning: WARNING: duplicate atom (ATOM     10  O   MET A   1      31.352   9.991 -20.514  1.00 14.01           O  , ATOM      3  O   MET A   1      33.432   9.978 -17.748  1.00  0.00           O  )
  warnings.warn("WARNING: duplicate atom (%s, %s)" % (atom, old_atom._pdb_string(old_atom.serial_number, atom.alternate_location_indicator)))
Adding hydrogens...
unmatched residues: [, ]
Traceback (most recent call last):
  File "00_simple.py", line 14, in 
    modeller.addHydrogens(omm_forcefield)
  File "/home/asgharrazavi/miniconda3/lib/python3.7/site-packages/simtk/openmm/app/modeller.py", line 943, in addHydrogens
    system = forcefield.createSystem(newTopology, rigidWater=False, nonbondedMethod=CutoffNonPeriodic)
  File "/home/asgharrazavi/miniconda3/lib/python3.7/site-packages/simtk/openmm/app/forcefield.py", line 1146, in createSystem
    raise ValueError('No template found for residue %d (%s).  %s' % (res.index+1, res.name, _findMatchErrors(self, res)))
ValueError: No template found for residue 1 (MET).  The set of atoms is similar to AOBT, but it is missing 2 atoms.

asgharrazavi · 2020-09-23T21:01:07Z

It appears that part of the problem might be in detecting hydrogen atoms. In the above PDB file, the first residue (capped MET) has 12 hydrogens, but the program detects 11.

# of hydrogens for residue : ['H', 'H2', 'H3', 'HA', 'HB2', 'HB3', 'HE1', 'HE2', 'HE3', 'HG2', 'HG3'] 11
# of hydrogens for residue : ['H', 'HA', 'HB', 'HG1', 'HG21', 'HG22', 'HG23'] 7
# of hydrogens for residue : ['H', 'HA', 'HB2', 'HB3', 'HG2', 'HG3'] 6
# of hydrogens for residue : ['H', 'H1', 'H2', 'H3'] 4

To get above information, I added the code print('# of hydrogens for residue %s:' %(residue), [i5.name for i5 in hydrogens], len(hydrogens)) to the modeller.py script:

                    hydrogens = [h for h in spec.hydrogens if (variant is None and pH <= h.maxph) or (h.variants is None and pH <= h.maxph) or (h.variants is not None and variant in h.variants)]
                    hydrogens = [h for h in hydrogens if h.terminal is None or (isNTerminal and h.terminal == 'N') or (isCTerminal and h.terminal == 'C')]
                    hydrogens = [h for h in hydrogens if h.parent in parentNames]

                    print('# of hydrogens for residue %s:' %(residue), [i5.name for i5 in hydrogens], len(hydrogens))   #ASGHAR 
                    # Loop over atoms in the residue, adding them to the new topology along with required hydrogens.

                    for parent in residue.atoms():

peastman · 2020-09-25T00:33:15Z

I think I see. Because the PDB reader expects ACE to be a separate residue, and doesn't add bonds correctly if it's merged. But CHARMM expects it to be part of the same residue, and doesn't work properly if it's separate.

Sigh.

Things would be simpler if everyone in the community would just agree on a consistent way of treating things like this. CHARMM is the nonstandard one in this case. The PDB lists ACE as a separate residue, and most programs follow it. But CHARMM calls it a patch instead of a residue, which requires it to be merged.

Anyone have suggestions on a clean way of dealing with this?

jchodera · 2020-09-25T02:20:49Z

Tagging in @j-wags, since we're working on the Open Force Field Topology extension to biopolymers right now.

Our approach will be more in line with perceiving biopolymer residues using industry-standard SMARTS matches, and then making this information available for traditional template-based schemes.

asgharrazavi · 2020-09-25T19:17:28Z

For the moment, is there a way to manually give information to topology builder for the caped N- and C-terminal residues? There is only 2 of them, so, it will help a lot just to have a manual adjustment before the problem is fundamentally fixed. Thanks.

peastman · 2020-09-25T21:15:34Z

You can add the missing bonds to the Topology by calling addBond() on it.

asgharrazavi · 2020-09-25T21:30:04Z

I created an N-capped MET residue template like this:

omm_forcefield = app.ForceField('charmm36.xml', 'charmm36/water.xml')
met = omm_forcefield._templates['MET']
metn = copy.deepcopy(met)
atoms = [atom for atom in met.atoms]
cay = copy.deepcopy(atoms[0])
cy = copy.deepcopy(atoms[0])
oy = copy.deepcopy(atoms[0])
hy1 = copy.deepcopy(atoms[0])
hy2 = copy.deepcopy(atoms[0])
hy3 = copy.deepcopy(atoms[0])

hy1.parameters, hy1.name, hy1.type, =  {'charge':0.09 } ,"HY1" ,"HA3", 
hy2.parameters, hy2.name, hy2.type, =  {'charge':0.09 } ,"HY2" ,"HA3", 
hy3.parameters, hy3.name, hy3.type, =  {'charge':0.09 } ,"HY3" ,"HA3", 
cy.parameters, cy.name, cy.type,    =  {'charge':0.51 } ,"CY"  ,"C", 
oy.parameters, oy.name, oy.type,    =  {'charge':-0.51} ,"OY"  ,"O", 

metn.addAtom(cay)  #index 17
metn.addAtom(hy1)  # 18
metn.addAtom(hy2)  # 19
metn.addAtom(hy3)  # 20
metn.addAtom(cy)   # 21
metn.addAtom(oy)   # 22

metn.addBond(17, 18)
metn.addBond(17, 19)
metn.addBond(17, 20)
metn.addBond(17, 21)
metn.addBond(21, 22)
metn.addBond(21, 0)
omm_forcefield._templates['METN'] = metn

I am pretty sure all the extra bonds are added now, but still get the "no template found" error :(

peastman · 2020-09-28T18:21:20Z

Did you add the missing bonds to the topology like I suggested?

asgharrazavi · 2020-09-28T21:12:23Z

Oh, I think I just fixed it! As the program outputs as "warnings" when reading PDB file, there are duplicate names for some heavy atoms in the last residue (GLU). I changed the last residue patching atom name like this and now it is working without error.

ATOM     38  N   GLU A   3      27.186  11.203 -18.842  1.00  8.87           N
ATOM     39  CA  GLU A   3      25.982  11.973 -18.580  1.00  7.86           C
ATOM     40  C   GLU A   3      25.238  11.354 -17.409  1.00 11.09           C
ATOM     41  O   GLU A   3      25.183  10.124 -17.284  1.00 10.54           O
ATOM     42  CB  GLU A   3      25.055  11.978 -19.791  1.00  9.04           C
ATOM     43  CG  GLU A   3      25.462  12.978 -20.839  1.00 12.16           C
ATOM     44  CD  GLU A   3      24.509  13.090 -22.010  1.00 19.31           C
ATOM     45  OE1 GLU A   3      24.886  13.774 -22.981  1.00 11.19           O
ATOM     46  OE2 GLU A   3      23.403  12.500 -21.976  1.00 12.96           O1-
ATOM     47  H   GLU A   3      27.134  10.454 -19.517  1.00  0.00           H
ATOM     48  HA  GLU A   3      26.256  12.998 -18.332  1.00  0.00           H
ATOM     49  HB3 GLU A   3      24.036  12.188 -19.466  1.00  0.00           H
ATOM     50  HB2 GLU A   3      25.030  10.981 -20.232  1.00  0.00           H
ATOM     51  HG3 GLU A   3      26.458  12.730 -21.207  1.00  0.00           H
ATOM     52  HG2 GLU A   3      25.581  13.958 -20.376  1.00  0.00           H
ATOM     53  NT  GLU A   3      24.622  12.170 -16.482  1.00  0.00           N
ATOM     54  CAT GLU A   3      23.907  11.584 -15.358  1.00  0.00           C
ATOM     55  HT  GLU A   3      24.667  13.173 -16.584  1.00  0.00           H
ATOM     56 1HA  GLU A   3      23.976  10.475 -15.428  1.00  0.00           H
ATOM     57 2HA  GLU A   3      24.361  11.929 -14.402  1.00  0.00           H
ATOM     58 3HA  GLU A   3      22.839  11.897 -15.386  1.00  0.00           H

So, the lesson is not to ignore "warnings" :).

EDIT:

I also needed to change hydrogen atoms in first residue like this:

ATOM      4 1HA  MET A   1      35.403   8.426 -17.715  1.00  0.00           H
ATOM      5 2HA  MET A   1      35.362   8.313 -19.530  1.00  0.00           H
ATOM      6 3HA  MET A   1      36.359   9.596 -18.727  1.00  0.00           H

Having trouble with un-hiding some previous comments, so here is the final PDB file and the code:

working PDB:

CRYST1   41.604  116.842   91.486  90.00  90.00  90.00 C 2 2 21      8
MODEL        1
ATOM      1  CY3 MET A   1      35.416   9.009 -18.663  1.00  0.00           C  
ATOM      2  CY  MET A   1      34.218   9.958 -18.693  1.00  0.00           C  
ATOM      3  OY  MET A   1      33.432   9.978 -17.748  1.00  0.00           O  
ATOM      4 1HA  MET A   1      35.403   8.426 -17.715  1.00  0.00           H  
ATOM      5 2HA  MET A   1      35.362   8.313 -19.530  1.00  0.00           H  
ATOM      6 3HA  MET A   1      36.359   9.596 -18.727  1.00  0.00           H  
ATOM      7  N   MET A   1      34.018  10.791 -19.775  1.00 20.33           N  
ATOM      8  CA  MET A   1      32.868  11.699 -19.798  1.00 18.20           C  
ATOM      9  C   MET A   1      31.568  10.891 -19.713  1.00 15.92           C  
ATOM     10  O   MET A   1      31.352   9.991 -20.514  1.00 14.01           O  
ATOM     11  CB  MET A   1      32.903  12.559 -21.061  1.00 21.42           C  
ATOM     12  CG  MET A   1      31.963  13.738 -21.031  1.00 25.49           C  
ATOM     13  SD  MET A   1      32.638  15.208 -21.861  1.00 31.29           S  
ATOM     14  CE  MET A   1      32.598  14.706 -23.443  1.00 28.17           C  
ATOM     15  HN  MET A   1      34.664  10.775 -20.551  1.00  0.00           H  
ATOM     16  HA  MET A   1      32.926  12.355 -18.930  1.00  0.00           H  
ATOM     17  HB3 MET A   1      32.676  11.938 -21.927  1.00  0.00           H  
ATOM     18  HB2 MET A   1      33.920  12.912 -21.229  1.00  0.00           H  
ATOM     19  HG3 MET A   1      31.724  13.984 -19.996  1.00  0.00           H  
ATOM     20  HG2 MET A   1      31.018  13.460 -21.499  1.00  0.00           H  
ATOM     21  HE1 MET A   1      32.982  15.498 -24.086  1.00  0.00           H  
ATOM     22  HE2 MET A   1      33.215  13.815 -23.558  1.00  0.00           H  
ATOM     23  HE3 MET A   1      31.571  14.475 -23.726  1.00  0.00           H  
ATOM     24  N   THR A   2      30.712  11.223 -18.741  1.00  9.98           N  
ATOM     25  CA  THR A   2      29.469  10.501 -18.503  1.00  8.85           C  
ATOM     26  C   THR A   2      28.342  11.463 -18.217  1.00 11.00           C  
ATOM     27  O   THR A   2      28.506  12.405 -17.424  1.00  9.24           O  
ATOM     28  CB  THR A   2      29.651   9.536 -17.306  1.00 15.81           C  
ATOM     29  OG1 THR A   2      30.789   8.703 -17.531  1.00 13.83           O  
ATOM     30  CG2 THR A   2      28.415   8.664 -17.034  1.00 15.49           C  
ATOM     31  H   THR A   2      30.915  12.004 -18.134  1.00  0.00           H  
ATOM     32  HA  THR A   2      29.222   9.920 -19.391  1.00  0.00           H  
ATOM     33  HB  THR A   2      29.843  10.135 -16.416  1.00  0.00           H  
ATOM     34  HG1 THR A   2      30.901   8.105 -16.788  1.00  0.00           H  
ATOM     35 HG21 THR A   2      28.610   8.012 -16.183  1.00  0.00           H  
ATOM     36 HG22 THR A   2      27.560   9.303 -16.813  1.00  0.00           H  
ATOM     37 HG23 THR A   2      28.197   8.058 -17.913  1.00  0.00           H  
ATOM     38  N   GLU A   3      27.186  11.203 -18.842  1.00  8.87           N  
ATOM     39  CA  GLU A   3      25.982  11.973 -18.580  1.00  7.86           C  
ATOM     40  C   GLU A   3      25.238  11.354 -17.409  1.00 11.09           C  
ATOM     41  O   GLU A   3      25.183  10.124 -17.284  1.00 10.54           O  
ATOM     42  CB  GLU A   3      25.055  11.978 -19.791  1.00  9.04           C  
ATOM     43  CG  GLU A   3      25.462  12.978 -20.839  1.00 12.16           C  
ATOM     44  CD  GLU A   3      24.509  13.090 -22.010  1.00 19.31           C  
ATOM     45  OE1 GLU A   3      24.886  13.774 -22.981  1.00 11.19           O  
ATOM     46  OE2 GLU A   3      23.403  12.500 -21.976  1.00 12.96           O1-
ATOM     47  H   GLU A   3      27.134  10.454 -19.517  1.00  0.00           H  
ATOM     48  HA  GLU A   3      26.256  12.998 -18.332  1.00  0.00           H  
ATOM     49  HB3 GLU A   3      24.036  12.188 -19.466  1.00  0.00           H  
ATOM     50  HB2 GLU A   3      25.030  10.981 -20.232  1.00  0.00           H  
ATOM     51  HG3 GLU A   3      26.458  12.730 -21.207  1.00  0.00           H  
ATOM     52  HG2 GLU A   3      25.581  13.958 -20.376  1.00  0.00           H  
ATOM     53  NT  GLU A   3      24.622  12.170 -16.482  1.00  0.00           N  
ATOM     54  CAT GLU A   3      23.907  11.584 -15.358  1.00  0.00           C  
ATOM     55  HT  GLU A   3      24.667  13.173 -16.584  1.00  0.00           H  
ATOM     56 1HA  GLU A   3      23.976  10.475 -15.428  1.00  0.00           H  
ATOM     57 2HA  GLU A   3      24.361  11.929 -14.402  1.00  0.00           H  
ATOM     58 3HA  GLU A   3      22.839  11.897 -15.386  1.00  0.00           H  
ENDMDL
END

Code:

from simtk.openmm.app import PDBFile
from simtk.openmm import app

print('building protein and solvating...')

# Load the force field parameters through OpenMM.
omm_forcefield = app.ForceField('charmm36.xml')

# Load the PDB file.
prot_pdb_file = PDBFile('test.pdb') 
modeller = app.Modeller(prot_pdb_file.topology, prot_pdb_file.positions)
top = prot_pdb_file.topology


## ADDING EXTRA BONDS ##
# Adding extra bonds to the first residue
print('Total number of bonds before adding extra bonds:', top.bonds)
res = [res for res in top.residues()][0]
atoms = [atom for atom in res.atoms()]
bonds = [bond for bond in res.bonds()]

print('Bonds for residue 1 before adding extra bonds to topology:')
for bond in bonds:
    print('%s %s' %(bond[0],bond[1]))

top.addBond(atoms[0],atoms[3])                                                                                                                                                                       
top.addBond(atoms[0],atoms[4])                                                                                                                                                                       
top.addBond(atoms[0],atoms[5])                                                                                                                                                                       
top.addBond(atoms[0],atoms[1])                                                                                                                                                                       
top.addBond(atoms[1],atoms[2])                                                                                                                                                                       
top.addBond(atoms[1],atoms[6])   
res = [res for res in top.residues()][0]
bonds = [bond for bond in res.bonds()]

print('Bonds for residue 1 after adding extra bonds to topology:')
for bond in bonds:
    print('%s %s' %(bond[0],bond[1]))


# Adding extra bonds to the last residue
print('Total number of bonds before adding extra bonds to last residue:', top.bonds)
res = [res for res in top.residues()][-1]
atoms = [atom for atom in top.atoms()]
bonds = [bond for bond in res.bonds()]

print('Bonds for the last residue before adding extra bonds to topology:')
for bond in bonds:
    print('%s %s' %(bond[0],bond[1]))

top.addBond(atoms[39],atoms[52])
top.addBond(atoms[52],atoms[53])
top.addBond(atoms[52],atoms[54])
top.addBond(atoms[53],atoms[55])
top.addBond(atoms[53],atoms[56])
top.addBond(atoms[53],atoms[57])
res = [res for res in top.residues()][-1]
bonds = [bond for bond in res.bonds()]

print('Bonds for the last residue after adding extra bonds to topology:')
for bond in bonds:
    print('%s %s' %(bond[0],bond[1]))
print('Total number of bonds after adding extra bonds:', top.bonds)
## FINISHED ADDING EXTRA BONDS ##


print('Adding hydrogens...')
modeller.addHydrogens(omm_forcefield)

XiaojuanHu · 2020-09-29T17:17:46Z

Hi, Thank you for the solution! I followed this method, but my code still didn't work. My system is Ac-Glu-NMe-Ca2+. Here is my pdb file:

HEADER    Dipeptide in gas
MODEL        0
ATOM      1 CAY  GLU     1       0.337  -0.068   0.143  1.00  0.00          C 
ATOM      2 CY   GLU     1      -0.082   0.185   1.565  1.00  0.00          C 
ATOM      3 OY   GLU     1       0.735   0.658   2.385  1.00  0.00          O 
ATOM      4 HY1  GLU     1      -0.461  -0.469  -0.490  1.00  0.00          H 
ATOM      5 HY2  GLU     1       1.177  -0.775   0.150  1.00  0.00          H 
ATOM      6 HY3  GLU     1       0.705   0.869  -0.293  1.00  0.00          H 
ATOM      7 N    GLU     1      -1.358  -0.135   1.883  1.00  0.00          N 
ATOM      8 CA   GLU     1      -2.015  -0.103   3.198  1.00  0.00          C 
ATOM      9 C    GLU     1      -2.171   1.347   3.703  1.00  0.00          C 
ATOM     10 O    GLU     1      -1.320   2.240   3.480  1.00  0.00          O 
ATOM     11 H    GLU     1      -1.916  -0.506   1.122  1.00  0.00          H 
ATOM     12 HA   GLU     1      -3.036  -0.450   2.992  1.00  0.00          H 
ATOM     13 CB   GLU     1      -1.374  -1.143   4.196  1.00  0.00          C 
ATOM     14 CG   GLU     1      -1.216  -0.767   5.687  1.00  0.00          C 
ATOM     15 CD   GLU     1      -0.148   0.309   5.796  1.00  0.00          C 
ATOM     16 OE1  GLU     1       1.019   0.035   5.348  1.00  0.00          O 
ATOM     17 OE2  GLU     1      -0.475   1.506   6.087  1.00  0.00          O 
ATOM     18 HB1  GLU     1      -0.372  -1.386   3.823  1.00  0.00          H 
ATOM     19 HB2  GLU     1      -1.973  -2.059   4.120  1.00  0.00          H 
ATOM     20 HG1  GLU     1      -0.879  -1.670   6.212  1.00  0.00          H 
ATOM     21 HG2  GLU     1      -2.149  -0.429   6.152  1.00  0.00          H 
ATOM     22 NT   GLU     1      -3.279   1.598   4.403  1.00  0.00          N 
ATOM     23 CAT  GLU     1      -3.511   2.868   5.086  1.00  0.00          C 
ATOM     24 HNT  GLU     1      -3.950   0.849   4.538  1.00  0.00          H 
ATOM     25 HT1  GLU     1      -3.436   3.698   4.375  1.00  0.00          H 
ATOM     26 HT2  GLU     1      -2.768   3.003   5.882  1.00  0.00          H 
ATOM     27 HT3  GLU     1      -4.515   2.853   5.517  1.00  0.00          H 
ATOM     28 CAL  CAL     2       0.822   1.959   4.229  1.00  0.00          Ca
END

My code is:

from simtk.openmm.app import *
from simtk.openmm.app.internal import *
from simtk.openmm import *
from simtk.unit import *
from sys import stdout
import os, sys, re

from xml.etree import ElementTree as ET

pdb = PDBFile('input_charmm.pdb')
forcefield = ForceField('charmm36.xml', 'charmm36/tip3p-pme-b.xml')

top = pdb.topology

## ADDING EXTRA BONDS ##
print('Total number of bonds before adding extra bonds:', top.bonds)
res = [res for res in top.residues()][0]
atoms = [atom for atom in res.atoms()]
bonds = [bond for bond in res.bonds()]

print('Bonds for residue 1 before adding extra bonds to topology:')
for bond in bonds:
	print('%s %s' %(bond[0],bond[1]))

top.addBond(atoms[0],atoms[3])                                                                                                                                                                       
top.addBond(atoms[0],atoms[4])                                                                                                                                                                       
top.addBond(atoms[0],atoms[5])                                                                                                                                                                       
top.addBond(atoms[0],atoms[1])                                                                                                                                                                       
top.addBond(atoms[1],atoms[2])                                                                                                                                                                       
top.addBond(atoms[1],atoms[6])
top.addBond(atoms[8],atoms[21])
top.addBond(atoms[21],atoms[22])
top.addBond(atoms[21],atoms[23])
top.addBond(atoms[22],atoms[24])
top.addBond(atoms[22],atoms[25])
top.addBond(atoms[22],atoms[26])   
res = [res for res in top.residues()][0]
bonds = [bond for bond in res.bonds()]

print('Bonds for residue 1 after adding extra bonds to topology:')
for bond in bonds:
    print('%s %s' %(bond[0],bond[1]))

## FINISHED ADDING EXTRA BONDS ##
system = forcefield.createSystem(top, nonbondedMethod=NoCutoff, constraints=None, removeCMMotion=False)

There are strange bonds appeared:
'Bonds for residue 1 before adding extra bonds to topology:
<Atom 8 (C) of chain 0 residue 0 (GLU)> <Atom 7 (CA) of chain 0 residue 0 (GLU)>
<Atom 8 (C) of chain 0 residue 0 (GLU)> <Atom 9 (O) of chain 0 residue 0 (GLU)>
<Atom 7 (CA) of chain 0 residue 0 (GLU)> <Atom 12 (CB) of chain 0 residue 0 (GLU)>
<Atom 7 (CA) of chain 0 residue 0 (GLU)> <Atom 11 (HA) of chain 0 residue 0 (GLU)>
<Atom 7 (CA) of chain 0 residue 0 (GLU)> <Atom 6 (N) of chain 0 residue 0 (GLU)>
<Atom 12 (CB) of chain 0 residue 0 (GLU)> <Atom 13 (CG) of chain 0 residue 0 (GLU)>
<Atom 12 (CB) of chain 0 residue 0 (GLU)> <Atom 18 (HB2) of chain 0 residue 0 (GLU)>
<Atom 12 (CB) of chain 0 residue 0 (GLU)> <Atom 17 (HB3) of chain 0 residue 0 (GLU)>
<Atom 14 (CD) of chain 0 residue 0 (GLU)> <Atom 13 (CG) of chain 0 residue 0 (GLU)>
<Atom 14 (CD) of chain 0 residue 0 (GLU)> <Atom 15 (OE1) of chain 0 residue 0 (GLU)>
<Atom 14 (CD) of chain 0 residue 0 (GLU)> <Atom 16 (OE2) of chain 0 residue 0 (GLU)>
<Atom 13 (CG) of chain 0 residue 0 (GLU)> <Atom 20 (HG2) of chain 0 residue 0 (GLU)>
<Atom 13 (CG) of chain 0 residue 0 (GLU)> <Atom 19 (HG3) of chain 0 residue 0 (GLU)>
<Atom 24 (H) of chain 0 residue 0 (GLU)> <Atom 6 (N) of chain 0 residue 0 (GLU)>
<Atom 25 (H2) of chain 0 residue 0 (GLU)> <Atom 6 (N) of chain 0 residue 0 (GLU)>
<Atom 26 (H3) of chain 0 residue 0 (GLU)> <Atom 6 (N) of chain 0 residue 0 (GLU)>'

Do I need to build a new template of capped Glu?

peastman · 2020-09-29T22:18:52Z

It recognizes the atom names "HT1", "HT2", and "HT3" as synonyms for "H", "H2", and "H3". And since those atoms would be bonded to N in a normal residue, it adds those bonds. Try giving them different names. @asgharrazavi called them "1HA", "2HA", and "3HA" in the example above.

XiaojuanHu · 2020-09-30T08:00:31Z

It worked! Thank you!

asgharrazavi · 2020-09-30T20:05:41Z

I think there is still an issue here. When I try to write the final PDB file using a suggestion from here, I see that there are extra hydrogen atoms:

code to write PDB:

with open('output.pdb', 'w') as f:
    PDBFile.writeFile(modeller.topology, modeller.positions, f)

output PDB for the first residue:

REMARK   1 CREATED WITH OPENMM 7.4.1, 2020-09-30
CRYST1   41.604  116.842   91.486  90.00  90.00  90.00 P 1           1 
ATOM      1  CY3 MET A   1      35.416   9.009 -18.663  1.00  0.00           C  
ATOM      2  CY  MET A   1      34.218   9.958 -18.693  1.00  0.00           C  
ATOM      3  OY  MET A   1      33.432   9.978 -17.748  1.00  0.00           O  
ATOM      4 1HY  MET A   1      35.394   8.440 -17.706  1.00  0.00           H  
ATOM      5 2HY  MET A   1      35.363   8.279 -19.498  1.00  0.00           H  
ATOM      6 3HY  MET A   1      36.373   9.568 -18.712  1.00  0.00           H  
ATOM      7  N   MET A   1      34.018  10.791 -19.775  1.00  0.00           N  
ATOM      8  H2  MET A   1      33.930  10.074 -20.521  1.00  0.00           H  
ATOM      9  H3  MET A   1      33.937  10.265 -18.847  1.00  0.00           H  
ATOM     10  CA  MET A   1      32.868  11.699 -19.798  1.00  0.00           C  
ATOM     11  C   MET A   1      31.568  10.891 -19.713  1.00  0.00           C  
ATOM     12  O   MET A   1      31.352   9.991 -20.514  1.00  0.00           O  
ATOM     13  CB  MET A   1      32.903  12.559 -21.061  1.00  0.00           C  
ATOM     14  CG  MET A   1      31.963  13.738 -21.031  1.00  0.00           C  
ATOM     15  SD  MET A   1      32.638  15.208 -21.861  1.00  0.00           S  
ATOM     16  CE  MET A   1      32.598  14.706 -23.443  1.00  0.00           C  
ATOM     17  H   MET A   1      34.933  11.267 -19.802  1.00  0.00           H  
ATOM     18  HA  MET A   1      32.921  12.342 -18.927  1.00  0.00           H  
ATOM     19  HB3 MET A   1      32.698  11.950 -21.974  1.00  0.00           H  
ATOM     20  HB2 MET A   1      33.943  12.956 -21.175  1.00  0.00           H  
ATOM     21  HG3 MET A   1      31.754  14.026 -19.976  1.00  0.00           H  
ATOM     22  HG2 MET A   1      30.979  13.460 -21.477  1.00  0.00           H  
ATOM     23  HE1 MET A   1      32.887  15.535 -24.124  1.00  0.00           H  
ATOM     24  HE2 MET A   1      33.297  13.863 -23.624  1.00  0.00           H  
ATOM     25  HE3 MET A   1      31.577  14.377 -23.733  1.00  0.00           H

These are the extra hydrogens:

ATOM      8  H2  MET A   1      33.930  10.074 -20.521  1.00  0.00           H  
ATOM      9  H3  MET A   1      33.937  10.265 -18.847  1.00  0.00           H

I tried renaming H atoms for N-terminal patching for different names, but wasn't successful. Strangely, the patched last residue is fine.

asgharrazavi · 2020-10-01T14:40:48Z

I deleted the extra atoms using delete method. However, my simulation is crashing after minimization.

PDB:

CRYST1   41.604  116.842   91.486  90.00  90.00  90.00 C 2 2 21      8
MODEL        1
ATOM      1  CY3 MET A   1      35.416   9.009 -18.663  1.00  0.00           C  
ATOM      2  CY  MET A   1      34.218   9.958 -18.693  1.00  0.00           C  
ATOM      3  OY  MET A   1      33.432   9.978 -17.748  1.00  0.00           O  
ATOM      4 1HY  MET A   1      35.403   8.426 -17.715  1.00  0.00           H  
ATOM      5 2HY  MET A   1      35.362   8.313 -19.530  1.00  0.00           H  
ATOM      6 3HY  MET A   1      36.359   9.596 -18.727  1.00  0.00           H  
ATOM      7  N   MET A   1      34.018  10.791 -19.775  1.00 20.33           N  
ATOM      8  CA  MET A   1      32.868  11.699 -19.798  1.00 18.20           C  
ATOM      9  C   MET A   1      31.568  10.891 -19.713  1.00 15.92           C  
ATOM     10  O   MET A   1      31.352   9.991 -20.514  1.00 14.01           O  
ATOM     11  CB  MET A   1      32.903  12.559 -21.061  1.00 21.42           C  
ATOM     12  CG  MET A   1      31.963  13.738 -21.031  1.00 25.49           C  
ATOM     13  SD  MET A   1      32.638  15.208 -21.861  1.00 31.29           S  
ATOM     14  CE  MET A   1      32.598  14.706 -23.443  1.00 28.17           C  
ATOM     15  HN  MET A   1      34.664  10.775 -20.551  1.00  0.00           H  
ATOM     16  HA  MET A   1      32.926  12.355 -18.930  1.00  0.00           H  
ATOM     17  HB3 MET A   1      32.676  11.938 -21.927  1.00  0.00           H  
ATOM     18  HB2 MET A   1      33.920  12.912 -21.229  1.00  0.00           H  
ATOM     19  HG3 MET A   1      31.724  13.984 -19.996  1.00  0.00           H  
ATOM     20  HG2 MET A   1      31.018  13.460 -21.499  1.00  0.00           H  
ATOM     21  HE1 MET A   1      32.982  15.498 -24.086  1.00  0.00           H  
ATOM     22  HE2 MET A   1      33.215  13.815 -23.558  1.00  0.00           H  
ATOM     23  HE3 MET A   1      31.571  14.475 -23.726  1.00  0.00           H  
ATOM     24  N   THR A   2      30.712  11.223 -18.741  1.00  9.98           N  
ATOM     25  CA  THR A   2      29.469  10.501 -18.503  1.00  8.85           C  
ATOM     26  C   THR A   2      28.342  11.463 -18.217  1.00 11.00           C  
ATOM     27  O   THR A   2      28.506  12.405 -17.424  1.00  9.24           O  
ATOM     28  CB  THR A   2      29.651   9.536 -17.306  1.00 15.81           C  
ATOM     29  OG1 THR A   2      30.789   8.703 -17.531  1.00 13.83           O  
ATOM     30  CG2 THR A   2      28.415   8.664 -17.034  1.00 15.49           C  
ATOM     31  H   THR A   2      30.915  12.004 -18.134  1.00  0.00           H  
ATOM     32  HA  THR A   2      29.222   9.920 -19.391  1.00  0.00           H  
ATOM     33  HB  THR A   2      29.843  10.135 -16.416  1.00  0.00           H  
ATOM     34  HG1 THR A   2      30.901   8.105 -16.788  1.00  0.00           H  
ATOM     35 HG21 THR A   2      28.610   8.012 -16.183  1.00  0.00           H  
ATOM     36 HG22 THR A   2      27.560   9.303 -16.813  1.00  0.00           H  
ATOM     37 HG23 THR A   2      28.197   8.058 -17.913  1.00  0.00           H  
ATOM     38  N   GLU A   3      27.186  11.203 -18.842  1.00  8.87           N  
ATOM     39  CA  GLU A   3      25.982  11.973 -18.580  1.00  7.86           C  
ATOM     40  C   GLU A   3      25.238  11.354 -17.409  1.00 11.09           C  
ATOM     41  O   GLU A   3      25.183  10.124 -17.284  1.00 10.54           O  
ATOM     42  CB  GLU A   3      25.055  11.978 -19.791  1.00  9.04           C  
ATOM     43  CG  GLU A   3      25.462  12.978 -20.839  1.00 12.16           C  
ATOM     44  CD  GLU A   3      24.509  13.090 -22.010  1.00 19.31           C  
ATOM     45  OE1 GLU A   3      24.886  13.774 -22.981  1.00 11.19           O  
ATOM     46  OE2 GLU A   3      23.403  12.500 -21.976  1.00 12.96           O1-
ATOM     47  H   GLU A   3      27.134  10.454 -19.517  1.00  0.00           H  
ATOM     48  HA  GLU A   3      26.256  12.998 -18.332  1.00  0.00           H  
ATOM     49  HB3 GLU A   3      24.036  12.188 -19.466  1.00  0.00           H  
ATOM     50  HB2 GLU A   3      25.030  10.981 -20.232  1.00  0.00           H  
ATOM     51  HG3 GLU A   3      26.458  12.730 -21.207  1.00  0.00           H  
ATOM     52  HG2 GLU A   3      25.581  13.958 -20.376  1.00  0.00           H  
ATOM     53  NT  GLU A   3      24.622  12.170 -16.482  1.00  0.00           N  
ATOM     54  CAT GLU A   3      23.907  11.584 -15.358  1.00  0.00           C  
ATOM     55  HT  GLU A   3      24.667  13.173 -16.584  1.00  0.00           H  
ATOM     56 1HA  GLU A   3      23.976  10.475 -15.428  1.00  0.00           H  
ATOM     57 2HA  GLU A   3      24.361  11.929 -14.402  1.00  0.00           H  
ATOM     58 3HA  GLU A   3      22.839  11.897 -15.386  1.00  0.00           H  
ENDMDL
END

Code:

from simtk.openmm.app import PDBFile
from simtk.openmm import app
from simtk import unit
from simtk.openmm import Platform
from simtk.openmm import LangevinIntegrator

print('building protein and solvating...')
# Load the force field parameters through OpenMM.
omm_forcefield = app.ForceField('charmm36.xml','charmm36/water.xml')

# Load the PDB file.
prot_pdb_file = PDBFile('test.pdb') 
modeller = app.Modeller(prot_pdb_file.topology, prot_pdb_file.positions)
top = prot_pdb_file.topology

## ADDING EXTRA BONDS ##
# Adding extra bonds to the first residue
res = [res for res in top.residues()][0]
atoms = [atom for atom in res.atoms()]
bonds = [bond for bond in res.bonds()]
top.addBond(atoms[0],atoms[3])                                                                                                                          
top.addBond(atoms[0],atoms[4])                                                                                                                          
top.addBond(atoms[0],atoms[5])                                                                                                                          
top.addBond(atoms[0],atoms[1])                                                                                                                          
top.addBond(atoms[1],atoms[2])                                                                                                                          
top.addBond(atoms[1],atoms[6])   

res = [res for res in top.residues()][0]
bonds = [bond for bond in res.bonds()]

# Adding extra bonds to the last residue
res = [res for res in top.residues()][-1]
atoms = [atom for atom in top.atoms()]
bonds = [bond for bond in res.bonds()]
top.addBond(atoms[39],atoms[52])
top.addBond(atoms[52],atoms[53])
top.addBond(atoms[52],atoms[54])
top.addBond(atoms[53],atoms[55])
top.addBond(atoms[53],atoms[56])
top.addBond(atoms[53],atoms[57])
res = [res for res in top.residues()][-1]
bonds = [bond for bond in res.bonds()]
## FINISHED ADDING EXTRA BONDS ##

print('Adding hydrogens...')
modeller.addHydrogens(omm_forcefield)

# Deleting exra hydrogens
h2 = [atom for atom in modeller.topology.atoms() if atom.name == 'H2' and atom.residue.name == 'MET']
print('atom to be deleted:', h2)
modeller.delete(h2)
h3 = [atom for atom in modeller.topology.atoms() if atom.name == 'H3' and atom.residue.name == 'MET']
print('atom to be deleted:', h3)
modeller.delete(h3)

print('Adding solvent...')
# here we can also add ionic concentration (read addSolvent?)
modeller.addSolvent(omm_forcefield, model='tip3p', padding=1*unit.nanometer, ionicStrength=0.1*unit.molar, positiveIon='Na+', negativeIon='Cl-')

with open('output.pdb', 'w') as f:
    PDBFile.writeFile(modeller.topology, modeller.positions, f)

prot_system = omm_forcefield.createSystem(modeller.topology, nonbondedMethod=app.PME, rigidWater=True)

## simulations
dt = 0.005*unit.picoseconds
temperature = 310*unit.kelvin
friction = 1.0/unit.picosecond
constraintTolerance = 0.000001
log_name = 'log.log'

# Simulation Options
steps = 160000
equilibrationSteps = 0
pdbReporter = app.PDBReporter('test.pdb', 20)
dataReporter = app.StateDataReporter(log_name, 1, totalSteps=steps, step=True, time=True, speed=True, progress=True, elapsedTime=True, remainingTime=True, potentialEnergy=True, kineticEnergy=True, totalEnergy=True, temperature=True, volume=True, density=True, separator=',')

integrator = LangevinIntegrator(temperature, friction, dt)
integrator.setConstraintTolerance(constraintTolerance)
simulation = app.Simulation(modeller.topology, prot_system, integrator)
simulation.context.setPositions(modeller.positions)
simulation.reporters.append(pdbReporter)
simulation.reporters.append(dataReporter)
simulation.minimizeEnergy()

print('writing minimized structure...')
positions = simulation.context.getState(getPositions=True).getPositions()
PDBFile.writeFile(simulation.topology, positions, open('minimized.pdb', 'w'))

# Set velocity and loadCheckpoint if available
simulation.context.setVelocitiesToTemperature(temperature)
simulation.currentStep = 0

# Simulate
print('Simulating...')
simulation.step(steps)

Output:

building protein and solvating...
Adding hydrogens...
before applying patches: [template, matches]: [None, None]
after applying patches: [template, matches]: [None, None]

atom to be deleted: [<Atom 7 (H2) of chain 0 residue 0 (MET)>]
atom to be deleted: [<Atom 7 (H3) of chain 0 residue 0 (MET)>]

Adding solvent...
writing minimized structure...
Simulating...
Traceback (most recent call last):
  File "final_code_on_git_edit2.py", line 96, in 
    simulation.step(steps)
  File "/home/asgharrazavi/miniconda3/lib/python3.7/site-packages/simtk/openmm/app/simulation.py", line 132, in step
    self._simulate(endStep=self.currentStep+steps)
  File "/home/asgharrazavi/miniconda3/lib/python3.7/site-packages/simtk/openmm/app/simulation.py", line 235, in _simulate
    self._generate_reports(unwrapped, False)
  File "/home/asgharrazavi/miniconda3/lib/python3.7/site-packages/simtk/openmm/app/simulation.py", line 254, in _generate_reports
    reporter.report(self, state)
  File "/home/asgharrazavi/miniconda3/lib/python3.7/site-packages/simtk/openmm/app/statedatareporter.py", line 195, in report
    self._checkForErrors(simulation, state)
  File "/home/asgharrazavi/miniconda3/lib/python3.7/site-packages/simtk/openmm/app/statedatareporter.py", line 347, in _checkForErrors
    raise ValueError('Energy is NaN')
ValueError: Energy is NaN

I wonder if I am adding solvent properly?

asgharrazavi · 2020-10-01T18:53:46Z

So, the problem was my MD time step. I was inappropriately using 5 fs for simulation time step without hydrogen mass repartitioning. It works with 2 fs time step :)

aniketmoha9 · 2021-08-03T07:53:00Z

I am trying to simulate a 2 amino acid sequence using amber force field in openmm. But there is always an error
'No template found for residue 1 (PHE). The set of atoms is similar to NTYR, but it is missing 1 hydrogen atoms.'

`from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
from sys import stdout

pdb = PDBFile('Leu-Phe.pdb')
forcefield = ForceField('amber14-all.xml', 'amber14/protein.ff14SB.xml')
system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME, nonbondedCutoff=1nanometer, constraints=HBonds)
integrator = VerletIntegrator(0.001picoseconds)
simulation = Simulation(pdb.topology, system, integrator)
simulation.context.setPositions(pdb.positions)
simulation.minimizeEnergy(maxIterations=100)
print('Saving...')
positions = simulation.context.getState(getPositions=True).getPositions()
PDBFile.writeFile(simulation.topology, positions, open('output.pdb', 'w'))
print('Done')`

i have also attached the pdb file.

Can you comment what is the problem and how to resolve it. @peastman

peastman · 2021-08-03T17:53:01Z

What order are they supposed to be in? The PHE is listed before the LEU, but the reside numbers are the other way around: PHE is residue 2 and LEU is residue 1. And the PHE is terminated (notice the OXT and HXT atoms), as if it were at the end of a chain, while the LEU is not terminated.

aniketmoha9 · 2021-08-04T04:49:12Z

I created new pdb file of only Phe residue and tried to simulate that using openmm but still there is error coming.

Here is the pdb file

i was able to minimize the same Phe residue using rdkit MM Force Field but in openmm with amber force field i am not able to get the minimize conformation. Please let me know if i am using the right xml file for my force field or is there a problem with my pdb file.

jchodera · 2021-08-04T04:53:41Z

Could you attach a ZIPped version of the PDB file to your issue?

From a quick inspection of the image of the PDB file contents, it looks like the N and C terminus have just been truncated---natural amino acids end in zwitterionic termini, with NH3(+)- at the N-terminus and -CO2(-) at the C-terminus. You only have one proton at the N-terminus and one oxygen at the C-terminus, suggesting you clipped out an internal PHE residue and left dangling bonds. OpenMM won't know how to match that.

aniketmoha9 · 2021-08-04T05:06:30Z

The first one is a normal Phe amino acid with the terminals, the second one we have removed the OH from the C terminal and H from the N terminal. Please let me know how to fix the pdb files.
Phe_with_terminals.zip
Phe_clipped_terminals.zip
i have also tried the zwitter ion form of the Phe but still there is an error.

Here is the zwiiter ion pdb file.
Phe_ion.zip

peastman · 2021-08-04T16:02:51Z

Try using OpenMM Setup. It can fix improperly terminated chains.

aniketmoha9 · 2021-08-04T16:39:32Z

Thanks for the reply. I have been also using the openmm setup but the same error is coming.

aniketmoha9 · 2021-08-06T05:14:23Z

Can you attach a peptide pdb file which can be simulated in openmm, Because i have been trying different amino acid sequence and every time there is some error.

peastman · 2021-08-06T15:46:20Z

Try this one: https://github.com/openmm/openmm/blob/master/examples/ala_ala_ala.pdb

aniketmoha9 · 2021-08-07T05:31:02Z

Thank you @peastman , i am able to simulate ala-ala-ala pdb file, i was getting error due to my pdb file, I want to know how to create the correct pdb file from a given amino acid sequence. Is there any script to do that. Because the way i am creating the pdb file i am not able to simulate it, but the connections are all right as i can see from its 3d visualization.
My file

your file

in your file every amino acid is starting with a nterminal and ending with cterminal which is quite logical, but in my case that is not the case, that is why the error is coming. Currently i am using rdkit to convert the smiles representation to mol file and then mol to pdb using open babel. If you can tell me how you are generating the pdb file it would be very helpful.

peastman · 2021-08-07T18:37:34Z

Your file is severely malformed. All the atoms for each residue need to go together. Your file begins by listing all the non-hydrogens for all residues, then follows with the hydrogens for all residues. RDKit just views your SMILES string as an arbitrary molecule, not realizing it's made up of residues whose atoms need to be kept together.

I don't know what tools people use for building peptides from scratch. It isn't something I've had to do much. Maybe someone else can suggest something.

aniketmoha9 · 2021-08-08T08:43:13Z

Thanks for replying. there is a function in rdkit MolFromSequence which generates correct mol format of the amino acid sequence which is only for peptides.

biosafetylvl5 closed this as completed Jun 19, 2018

This comment has been minimized.

Sign in to view

peastman mentioned this issue Sep 26, 2020

No template found for residue #2858

Closed

This comment has been minimized.

Sign in to view

asgharrazavi mentioned this issue Oct 2, 2020

Missing bond types when converting a system to GROMACS top file openforcefield/openff-toolkit#731

Closed

DanielWicz mentioned this issue Aug 13, 2021

ValueError: No template found for residue 37 (TYR). The set of atoms is similar to TTHY, but it is missing 2 atoms. #3108

Closed

ixy90 mentioned this issue Apr 24, 2022

no template found for residue/ coordinate "..." could not be represented in a width-8 field #3582

Closed

No template found for residue #2103

No template found for residue #2103

Comments

biosafetylvl5 commented Jun 18, 2018

biosafetylvl5 commented Jun 18, 2018

peastman commented Jun 19, 2018

biosafetylvl5 commented Jun 19, 2018

asgharrazavi commented Sep 21, 2020

peastman commented Sep 22, 2020

This comment has been minimized.

asgharrazavi commented Sep 22, 2020 • edited

peastman commented Sep 22, 2020

asgharrazavi commented Sep 22, 2020

asgharrazavi commented Sep 23, 2020

peastman commented Sep 25, 2020

jchodera commented Sep 25, 2020

asgharrazavi commented Sep 25, 2020

peastman commented Sep 25, 2020

asgharrazavi commented Sep 25, 2020

peastman commented Sep 28, 2020

This comment has been minimized.

This comment has been minimized.

asgharrazavi commented Sep 28, 2020 • edited

Having trouble with un-hiding some previous comments, so here is the final PDB file and the code:

working PDB:

Code:

XiaojuanHu commented Sep 29, 2020

peastman commented Sep 29, 2020

XiaojuanHu commented Sep 30, 2020

asgharrazavi commented Sep 30, 2020

asgharrazavi commented Oct 1, 2020 • edited

asgharrazavi commented Oct 1, 2020 • edited

aniketmoha9 commented Aug 3, 2021 • edited

peastman commented Aug 3, 2021

aniketmoha9 commented Aug 4, 2021

jchodera commented Aug 4, 2021

aniketmoha9 commented Aug 4, 2021 • edited

peastman commented Aug 4, 2021

aniketmoha9 commented Aug 4, 2021

aniketmoha9 commented Aug 6, 2021

peastman commented Aug 6, 2021

aniketmoha9 commented Aug 7, 2021 • edited

peastman commented Aug 7, 2021

aniketmoha9 commented Aug 8, 2021

asgharrazavi commented Sep 22, 2020 •

edited

asgharrazavi commented Sep 28, 2020 •

edited

asgharrazavi commented Oct 1, 2020 •

edited

asgharrazavi commented Oct 1, 2020 •

edited

aniketmoha9 commented Aug 3, 2021 •

edited

aniketmoha9 commented Aug 4, 2021 •

edited

aniketmoha9 commented Aug 7, 2021 •

edited