pdb2gmx generates wrong virtual sites for tetrahedral amine
Summary
I am to simulate a protein with a modified lysine, where the side chain has been extended so that the amine is now a secondary amine and neutral. It should be in a tetrahedral configuration, with one H and one lone pair (the latter not explicitly modelled at present). I have created an rtp entry for the residue and added all required new parameters to the forcefield files, but pdb2gmx makes the wrong virtual sites parameters for the amine hydrogen. Specifically, it generates a virtual_sites3 type 2 (3fd) entry, which makes the amine flat, but I would expect a virtual_site3 type 4 (3out). Changing it manually is error prone and tedious. Also this error could go unnoticed, as it did during my first round of simulations.
GROMACS version
gmx -quiet --version
GROMACS version: 2022.2
Precision: mixed
Memory model: 64 bit
MPI library: thread_mpi
OpenMP support: enabled (GMX_OPENMP_MAX_THREADS = 128)
GPU support: disabled
SIMD instructions: ARM_NEON_ASIMD
CPU FFT library: fftw-3.3.8
GPU FFT library: none
TNG support: enabled
Hwloc support: disabled
Tracing support: disabled
C compiler: /opt/homebrew/bin/gcc-11 GNU 11.2.0
C compiler flags: -Wno-missing-field-initializers -fexcess-precision=fast -funroll-all-loops -O3 -DNDEBUG
C++ compiler: /opt/homebrew/bin/g++-11 GNU 11.2.0
C++ compiler flags: -Wno-missing-field-initializers -fexcess-precision=fast -funroll-all-loops -fopenmp -O3 -DNDEBUG
Steps to reproduce
gmx pdb2gmx -ignh -water tip3p -f LFR_A.pdb -o LFR_p2g -i LFR_p2g -p LFR_p2g -vsite hydrogens
Choose amber99sb-ildn LFR/S.
I have tested using the resulting topology for an energy minimization (after solvating etc) and the amine becomes planar during the run. If I manually edit the topology (like the expected result below) I get a tetrahedral geometry. (Please describe how we can reproduce the bug, and share all files needed - ideally both the TPR file and the raw GRO/MDP/TOP files needed to regenerate it. Bugs that only appear after running for 3 hours on 200 GPUs unfortunately tend to not get a lot of attention. You will typically get much faster attention if you have been able to narrow it down to the smallest possible input, command line, system size, etc.)
What is the current bug behavior?
pdb2gmx makes virtual site parameters that makes a tetrahedral amine flat: a 3fd site.
[ virtual_sites3 ]
...
1450 1449 1446 1451 2
(1450 is the H on the amine)
What did you expect the correct behavior to be?
A tetrahedral geometry with 3out construction.
[ virtual_sites3 ]
...
1450 1449 1446 1451 4
Possible fixes
I have not yet found where the vsites are handled in pdb2gmx, but I can probably have a closer look if someone gives me a few pointers (no pun intended).