Polarizable CHARMM drude ff - missing parameters?

GROMACS version: 2016-dev-20181220-a2cab74156
GROMACS modification: Yes

Hello,

I am trying to setup a simulation using the drude-2013f_2018_10.ff available at: http://mackerell.umaryland.edu/charmm_drude_ff.shtml

I succesfully branched and installed the Gromacs version supporting polarisability (following the instructions at that webside) and successfully (or so it seems) created a topology using the modified version of pdb2gmx within that Gromacs installation.
The protein I’m trying to setup is the textbook example lysozime (1AKI).
I also explicitly set the termini residues to NTES (for charged N-terminus) and CTES (charged C-terminus).

My PDB2GMX command is: gmx_drude pdb2gmx -f 1AKI.pdb -p -o pdb2gmx_out -ter

PDB2GMX does its job, seemingly correctly. Adding the right atoms according to aminoacids.n.tdb and aminoacids.c.tdb.
However, the GROMPP step fails with a bunch of errors relating to missing (backbone) parameters from bond types to U-bond types.

Here an extract of the errors from GROMPP:

ERROR 1 [file topol.top, line 6479]:
No default Bond types
[…]
ERROR 6 [file topol.top, line 28416]:
No default U-B types
[…]
ERROR 9 [file topol.top, line 33591]:
No default Proper Dih. types
[…]

In total, for 1AKI, I get 14 errors.

I am wondering that this has something to do with the initial topology setup by PDB2GMX? Couldn’t find anything similar on the forum and was wondering if anybody else has experienced the same problem?

Thanks in advance for your time in answering my question! Any help will be enormously appreciated!

Cheers,
Davide

This is incorrect. NTES and CTES are only for zwitterionic, single amino acids. You are likely getting the wrong patching/atom type assignment. What you want are NTER and CTER.

Thank Justin. Very much appreciated.

NTER and CTER don’t show up as options. That’s what shows up for N-terminus:
Select start terminus type for LYS-1
0: NH3+
1: NTES
2: GNTS
3: None

and C-terminus:
Select end terminus type for ARG-128
0: COO-
1: COOH
2: CTES
3: CNES
4: None

I am guessing that NTER and CTER refer to choices NH3+ and COO- as charged termini of a chain. However, selecting those I get the following error in the creation of the topology:

Fatal error:
Could not find Drude when adding 1-3 Thole: CA-OT2

Thank you.
D.

Yes, that’s correct. Sorry, I spend most of my time in CHARMM these days :)

This shouldn’t happen, and I recall having fixed the code long ago to address this. Can you please provide:

  1. The full screen output of pdb2gmx, including the full error (with code line and commit hash - this will help me trace it)
  2. The compiler (and version) you used to build GROMACS

Hello Justin,

Many thanks for all the help with this:

This is the full screen from pdb2gmx:

gmx_drude pdb2gmx -f 1AKI_noter.pdb -p -o pdb2gmx_out -ter

Select the Force Field:
From current directory:
1: Drude-2013b polarizable force field
From ‘/Users/mercadantehome/progs/compiled/Gromacs2020-drude/share/gromacs/top’:
2: AMBER03 protein, nucleic AMBER94 (Duan et al., J. Comp. Chem. 24, 1999-2012, 2003)
3: AMBER94 force field (Cornell et al., JACS 117, 5179-5197, 1995)
4: AMBER96 protein, nucleic AMBER94 (Kollman et al., Acc. Chem. Res. 29, 461-469, 1996)
5: AMBER99 protein, nucleic AMBER94 (Wang et al., J. Comp. Chem. 21, 1049-1074, 2000)
6: AMBER99SB protein, nucleic AMBER94 (Hornak et al., Proteins 65, 712-725, 2006)
7: AMBER99SB-ILDN protein, nucleic AMBER94 (Lindorff-Larsen et al., Proteins 78, 1950-58, 2010)
8: AMBERGS force field (Garcia & Sanbonmatsu, PNAS 99, 2782-2787, 2002)
9: CHARMM27 all-atom force field (CHARM22 plus CMAP for proteins)
10: GROMOS96 43a1 force field
11: GROMOS96 43a2 force field (improved alkane dihedrals)
12: GROMOS96 45a3 force field (Schuler JCC 2001 22 1205)
13: GROMOS96 53a5 force field (JCC 2004 vol 25 pag 1656)
14: GROMOS96 53a6 force field (JCC 2004 vol 25 pag 1656)
15: GROMOS96 54a7 force field (Eur. Biophys. J. (2011), 40, 843-856, DOI: 10.1007/s00249-011-0700-9)
16: OPLS-AA/L all-atom force field (2001 aminoacid dihedrals)
1

Using the Drude-2013f_2018_10 force field in directory ./drude-2013f_2018_10.ff

Opening force field file ./drude-2013f_2018_10.ff/watermodels.dat

Select the Water Model:
1: SWM4-NDP 4-site polarizable water model with negatively charged Drude
2: SWM6 6-site polarizable water model with negatively charged Drude
3: None
1
Opening force field file ./drude-2013f_2018_10.ff/aminoacids.r2b
Opening force field file ./drude-2013f_2018_10.ff/dna.r2b
Opening force field file ./drude-2013f_2018_10.ff/rna.r2b
Warning: Residue ‘DA’ already present with type ‘DNA’ in database, ignoring new type ‘RNA’.Warning: Residue ‘DA’ already present with type ‘DNA’ in database, ignoring new type ‘RNA’.Warning: Residue ‘DA’ already present with type ‘DNA’ in database, ignoring new type ‘RNA’.Warning: Residue ‘DA’ already present with type ‘DNA’ in database, ignoring new type ‘RNA’.Warning: Residue ‘DC’ already present with type ‘DNA’ in database, ignoring new type ‘RNA’.Warning: Residue ‘DC’ already present with type ‘DNA’ in database, ignoring new type ‘RNA’.Warning: Residue ‘DC’ already present with type ‘DNA’ in database, ignoring new type ‘RNA’.Warning: Residue ‘DC’ already present with type ‘DNA’ in database, ignoring new type ‘RNA’.Warning: Residue ‘DG’ already present with type ‘DNA’ in database, ignoring new type ‘RNA’.Warning: Residue ‘DG’ already present with type ‘DNA’ in database, ignoring new type ‘RNA’.Warning: Residue ‘DG’ already present with type ‘DNA’ in database, ignoring new type ‘RNA’.Warning: Residue ‘DG’ already present with type ‘DNA’ in database, ignoring new type ‘RNA’.Reading 1AKI_noter.pdb…
WARNING: all CONECT records are ignored
Read 992 atoms
Analyzing pdb file
Splitting chemical chains based on TER records or chain id changing.
There are 1 chains and 0 blocks of water and 128 residues with 992 atoms

chain #res #atoms
1 ‘A’ 128 992

All occupancies are one
Opening force field file ./drude-2013f_2018_10.ff/atomtypes.atp
Atomtype 169
Reading residue database… (drude-2013f_2018_10)
Opening force field file ./drude-2013f_2018_10.ff/aminoacids.rtp
Detected a polarizable force field
Residue 253
Sorting it all out…
Opening force field file ./drude-2013f_2018_10.ff/dna.rtp
Detected a polarizable force field
Residue 261
Sorting it all out…
Opening force field file ./drude-2013f_2018_10.ff/rna.rtp
Detected a polarizable force field
Residue 266
Sorting it all out…
Opening force field file ./drude-2013f_2018_10.ff/aminoacids.hdb
Opening force field file ./drude-2013f_2018_10.ff/dna.hdb
Opening force field file ./drude-2013f_2018_10.ff/rna.hdb
Opening force field file ./drude-2013f_2018_10.ff/aminoacids.n.tdb
Opening force field file ./drude-2013f_2018_10.ff/dna.n.tdb
Opening force field file ./drude-2013f_2018_10.ff/rna.n.tdb
Opening force field file ./drude-2013f_2018_10.ff/aminoacids.c.tdb
Opening force field file ./drude-2013f_2018_10.ff/dna.c.tdb
Opening force field file ./drude-2013f_2018_10.ff/rna.c.tdb

Back Off! I just backed up topol.top to ./#topol.top.1#
Processing chain 1 ‘A’ (992 atoms, 128 residues)
Analysing hydrogen-bonding network for automated assignment of histidine
protonation. 212 donors and 183 acceptors were found.
There are 254 hydrogen bonds
Will use HISE for residue 15
Identified residue LYS1 as a starting terminus.
Identified residue ARG128 as a ending terminus.
8 out of 8 lines of specbond.dat converted successfully
Special Atom Distance matrix:
CYS6 MET12 HIS15 CYS30 CYS64 CYS76 CYS80
SG48 SD87 NE2118 SG238 SG513 SG601 SG630
MET12 SD87 1.166
HIS15 NE2118 1.776 1.019
CYS30 SG238 1.406 1.054 2.069
CYS64 SG513 2.835 1.794 1.789 2.241
CYS76 SG601 2.704 1.551 1.468 2.116 0.765
CYS80 SG630 2.959 1.951 1.916 2.391 0.199 0.944
CYS94 SG724 2.550 1.407 1.382 1.975 0.665 0.202 0.855
MET105 SD799 1.827 0.911 1.683 0.888 1.849 1.461 2.036
CYS115 SG889 1.576 1.084 2.078 0.200 2.111 1.989 2.262
CYS127 SG981 0.197 1.072 1.721 1.313 2.799 2.622 2.934
CYS94 MET105 CYS115
SG724 SD799 SG889
MET105 SD799 1.381
CYS115 SG889 1.853 0.790
CYS127 SG981 2.475 1.686 1.483
Linking CYS-6 SG-48 and CYS-127 SG-981…
Linking CYS-30 SG-238 and CYS-115 SG-889…
Linking CYS-64 SG-513 and CYS-80 SG-630…
Linking CYS-76 SG-601 and CYS-94 SG-724…
Select start terminus type for LYS-1
0: NH3+
1: NTES
2: GNTS
3: None
0
Start terminus LYS-1: NH3+
Select end terminus type for ARG-128
0: COO-
1: COOH
2: CTES
3: CNES
4: None
0
End terminus ARG-128: COO-
Opening force field file ./drude-2013f_2018_10.ff/aminoacids.arn
Opening force field file ./drude-2013f_2018_10.ff/dna.arn
Opening force field file ./drude-2013f_2018_10.ff/rna.arn
Checking for duplicate atoms…
Generating any missing hydrogen atoms, Drudes, lone pairs, and/or adding termini.
Now there are 128 residues with 3323 atoms
Making bonds…
Number of bonds was 2962, now 2959
Generating angles, dihedrals and pairs…
Generating Drude exclusions, lone pairs, and Thole screening…
Generating anisotropic polarization…wrote 195 entries.
Generating isotropic polarization…wrote 0 entries.
Generating virtual sites from lone pairs…wrote 373 virtual sites.


Program: gmx pdb2gmx, version 2016-dev-20181220-a2cab74156
Source file: src/gromacs/gmxpreprocess/gen_ad.cpp (line 1735)

Fatal error:
Could not find Drude when adding 1-3 Thole: CA-OT2

For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors

and this are the info on compilation (quick and “dirty” in this case):

GROMACS version: 2016-dev-20181220-a2cab74156
GIT SHA1 hash: a2cab7415652df7ca3631b5a750dddeda87573f9
Precision: single
Memory model: 64 bit
MPI library: MPI
OpenMP support: disabled
GPU support: disabled
SIMD instructions: AVX_256
FFT library: fftw-3.3.4-sse2
RDTSCP usage: enabled
TNG support: enabled
Tracing support: disabled
Built on: Fri 5 Jun 2020 21:26:09 NZST
Built by: mercadantehome@Davides-MBP [CMAKE]
Build OS/arch: Darwin 18.7.0 x86_64
Build CPU vendor: Unknown, detect failed
Build CPU brand: Unknown, detect failed
Build CPU family: 0 Model: 0 Stepping: 0
Build CPU features:
C compiler: /usr/local/bin/mpicc Clang 10.0.1.10010046
C compiler flags: -mavx -Wall -Wno-unused -Wunused-value -Wunused-parameter -Wno-unknown-pragmas -O3 -DNDEBUG
C++ compiler: /usr/local/bin/mpic++ Clang 10.0.1.10010046
C++ compiler flags: -mavx -std=c++0x -Wextra -Wno-missing-field-initializers -Wpointer-arith -Wall -Wno-unused-function -Wno-unknown-pragmas -O3 -DNDEBUG

Thanks. Looking forward to your reply.

Cheers,
Davide

I admit I’m stumped. In the past, I found that old (pre-5) versions of gcc generated incorrect code in this section though I never really resolved why. I’ve never tried compiling with the compilers you’ve used (and there’s no point assembling the entire gmx binary with MPI, anyway). There must be something upstream of this in the terminus patching that goes wrong. But admittedly all of the Drude construction on termini is a hack, because pdb2gmx doesn’t support adding atoms onto added atoms (e.g. it shouldn’t be possible to add DOT2 on top of OT2, which is added). I have a workaround in place, but this is something that is likely going to break pretty easily. It all works for me, though (Linux with gcc), so I thought the problem was fixed. Clearly it’s not.

Hello Justing,

thanks once again.
I have another machine (Ubuntu 18.04) and I have compiled the same code using GNU compilers 8.0.4 as I seemed to understand that compiling on Linux with gcc would solve the problem.
However, I still get the error.

What compiler version you are using?

Overall, is this going to be fixed at some stage or I should drop GMX altogether for calculations using a polarisation setup?

It seems that since GMX2019, polarisation (isotropic or anisotropic) is supported, however I am a bit confused how this relates to the hack in the drude branch.
http://manual.gromacs.org/documentation/current/reference-manual/functions/polarization.html

How does the topology get generated for applying polarisation in GMX 2019/2020?

Thanks for all the help.
Just still trying to understand if GMX is a good option for carrying out MD using polarisable ff.

Thank you.
D.

P.S. Please see below my compilation details:

GROMACS version: 2016-dev-20181220-a2cab74156-dirty
GIT SHA1 hash: a2cab7415652df7ca3631b5a750dddeda87573f9 (dirty)
Precision: single
Memory model: 64 bit
MPI library: thread_mpi
OpenMP support: enabled (GMX_OPENMP_MAX_THREADS = 32)
GPU support: disabled
SIMD instructions: AVX_256
FFT library: fftw-3.3.4-sse2-avx
RDTSCP usage: enabled
TNG support: enabled
Tracing support: disabled
Built on: Sun Jun 7 19:51:01 NZST 2020
Built by: davide@sc383507 [CMAKE]
Build OS/arch: Linux 5.3.0-53-generic x86_64
Build CPU vendor: AMD
Build CPU brand: AMD Ryzen 9 3900X 12-Core Processor
Build CPU family: 23 Model: 113 Stepping: 0
Build CPU features: aes apic avx avx2 clfsh cmov cx8 cx16 f16c fma htt lahf misalignsse mmx msr nonstop_tsc pclmuldq pdpe1gb popcnt pse rdrnd rdtscp sha sse2 sse3 sse4a sse4.1 sse4.2 ssse3
C compiler: /usr/bin/gcc-8 GNU 8.4.0
C compiler flags: -mavx -Wundef -Wextra -Wno-missing-field-initializers -Wno-sign-compare -Wpointer-arith -Wall -Wno-unused -Wunused-value -Wunused-parameter -O3 -DNDEBUG -funroll-all-loops -fexcess-precision=fast -Wno-array-bounds
C++ compiler: /usr/bin/g+±8 GNU 8.4.0
C++ compiler flags: -mavx -std=c++0x -Wundef -Wextra -Wno-missing-field-initializers -Wpointer-arith -Wall -Wno-unused-function -O3 -DNDEBUG -funroll-all-loops -fexcess-precision=fast -Wno-array-bounds

I am currently the only person working on this and my available time for coding is unfortunately very little. Now that I know something is broken that I thought I fixed about 7 months ago, I’ll probably go back into it more quickly. But I can’t promise when it will be ready. I’m aiming to include the Drude features in the 2021 release (which means I need to polish off the remaining code this summer, at least) but that’s a goal, not a promise. The core features all work (in my hands) in the Drude branch but you should treat anything from it as being beta testing. Don’t use it for science you hope to publish.

No relation at all. The current shell model support has largely been used for simulations of water. It doesn’t support everything that our more developed force field needs.

At present, I would say no, at least from the perspective of our Drude model. CHARMM, NAMD, and OpenMM are complete implementations. I apologize for the delay in getting our Drude stuff included in GROMACS (I know people want it) but again it’s just me doing the coding and only in spare time, which has been precious little for the last couple years.

Thanks Justin.

Got it, perfectly clear. While I hopefully wait to return to Gromacs once it will fully support polarisability via the Drude framework, I might (unwillingly) divert on other software (like OpenMM or NAMD for instance).

Of course no need to apologise, it would be great to see this implemented soon, but it’s totally understandable when you are the only one developing this!!

Thanks a lot for all the very appreciated answers!!

Cheers,
Davide