[gmx-users] Simulating an alkanethiol SAM on gold: atom type not recognized when running pdb2gmx

[Justin A. Lemkul]

Olivia Waring wrote:
> Hello GROMACS users,
> 
> I wish to simulate an alkanethiol self-assembled monolayer on a gold 
> surface. The first molecule I am attempting to model is quite simple: 
> CH3(CH2)11S. I have constructed a united atom force field based on 
> parameters published by Hautman and Klein (J. Chem. Phys. 1989) to 
> represent the alkanethiol chain. I wrote an rtp file containing a single 
> residue [ ALK ] which contains the following list of atoms:
> 
> [ ALK ]
>  [ atoms ]
>  ; name type charge   chargegroup
>    CH3  CH3   0.000   1
>    CH2A CH2   0.000   1
>    CH2B CH2   0.000   1
>    CH2C CH2   0.000   1
>    CH2D CH2   0.000   1
>    CH2E CH2   0.000   1
>    CH2F CH2   0.000   1
>    CH2G CH2   0.000   1
>    CH2H CH2   0.000   1
>    CH2I CH2   0.000   1
>    CH2J CH2   0.000   1
>    CH2K CH2   0.000   1
>    S    S    0.000   1
> 
> I have defined the CH3, CH2, and S atom types in an atomtypes.atp file. 
> I used PRODRG to generate an all-atom pdb file representing my 
> alkanethiol, which I then manually modified to reflect my united atom 
> approach. Here is an excerpt:
> 
> HETATM    1  CH3        ALK     1      -0.920  -1.820  -0.970  1.00 
> 20.00             C
> HETATM    2  CH2A       ALK     1      -0.010  -2.270   0.170  1.00 
> 20.00             C
> HETATM    3  CH2B       ALK     1       0.240  -3.780   0.060  1.00 
> 20.00             C
> HETATM    4  CH2C       ALK     1       1.140  -4.290   1.190  1.00 
> 20.00             C
> HETATM    5  CH2D       ALK     1       0.420  -4.310   2.540  1.00 
> 20.00             C
> HETATM    6  CH2E       ALK     1       1.310  -4.850   3.660  1.00 
> 20.00             C
> HETATM    7  CH2F       ALK     1       1.530  -6.360   3.540  1.00 
> 20.00             C
> HETATM    8  CH2G       ALK     1       2.420  -6.850   4.680  1.00 
> 20.00             C
> HETATM    9  CH2H       ALK     1       2.630  -8.360   4.540  1.00 
> 20.00             C
> HETATM   10  CH2I       ALK     1       3.510  -8.860   5.700  1.00 
> 20.00             C
> HETATM   11  CH2J       ALK     1       3.750 -10.370   5.610  1.00 
> 20.00             C
> HETATM   12  CH2K       ALK     1       2.490 -11.170   5.960  1.00 
> 20.00             C
> HETATM   12  S          ALK     1       2.730 -12.980   5.840  1.00 
> 20.00             S
> 
> When I try to run pdb2gmx on my coordinate file, I get the following error: 
> 
> Atom CH2 in residue ALK 0 was not found in rtp entry ALK with 13 
> atoms while sorting atoms.
> 

The error suggests something is wrong with the format of your .pdb file.  Likely 
a column is shifted incorrectly and the atom name is interpreted as CH2 rather 
than CH2[whatever].  Note that pdb2gmx is unnecessary in this case.  PRODRG 
provides a topology (albeit usually wrong) that can simply be #included in a 
system topology.  Generally, the charges are all that need to be adjusted, but 
if you're content with setting everything to zero, then the PRODRG topology is 
easy to fix.

> Although I have read through a significant portion of the manual, I am a 
> very new GROMACS user, and have most likely made numerous mistakes!! If 
> someone could please point out how to address this error and anything 
> else I might be missing, I would appreciate it very much. Any general 
> tips for constructing a system from scratch are also quite welcome.
> 

In a generic sense, it is as easy as:

#include "whatever.ff/forcefield.itp"

#include "molecule.itp"

#include "another_molecule.itp"

[ system ]
My system

[ molecules ]
MOL   1
MOL2  1

-Justin

-- 
========================================

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================