On Tue, May 22, 2018 at 7:57 AM, David A Case <david.case.rutgers.edu>
wrote:
> Hi everyone:
>
> This problem has come up repeatedly on the mailing list, and we have a
> simple test case, but I don't see a simple solution.
>
> Run the following with the attached files:
>
> pmemd.cuda -i md.in -p cA2.prmtop -c restrt
>
> You should get an inconsistent SHAKEH error, arising in gpu.cpp in the
> gpu_shake_setup() routine. The code seems to think that atom 62 is a
> hydrogen, but it is not.
>
> Some notes:
>
> 1. The code might be counting atoms from 0, rather than 1, but that
> would just make the reported error message misleading to users expecting
> (as parmed does) to have atom numbers starting at 1. But there aren't
> multiple bonds to hydrogens even if one interprets the atom numbers
> reported as being 0-based.
>
> 2. This only seems to happen with prmtops created by parmed/chamber.
> But the details of the atoms involved, including their masses, all look
> OK.
>
> 3. If you figure it out, please add some comments to this routine to
> make it easier to debug if a similar problem arises again.
>
This prmtop looks really weird. When I do some inspection with ParmEd, I
see the following:
> printDetails :4
The mask :4 matches 17 atoms:
ATOM RES RESNAME NAME TYPE At.# LJ Radius LJ Depth
Mass Charge GB Radius GB Screen
49 4 HSE N NH1 7 1.8500 0.2000
14.0070 -0.4700 1.5500 0.7900
50 4 HSE HN H 1 0.2245 0.0460
1.0080 0.3100 1.3000 0.8500
51 4 HSE CA CT1 6 2.2750 0.0200
12.0110 0.0700 1.7000 0.7200
52 4 HSE HA HB 1 1.3200 0.0220
1.0080 0.0900 1.3000 0.8500
53 4 HSE CB CT2 6 2.1750 0.0550
12.0110 -0.0800 1.7000 0.7200
54 4 HSE HB1 HA 1 1.3200 0.0220
1.0080 0.0900 1.3000 0.8500
55 4 HSE HB2 HA 1 1.3200 0.0220
1.0080 0.0900 1.3000 0.8500
56 4 HSE ND1 NR2 1 1.8500 0.2000
14.0070 -0.7000 1.3000 0.8500
57 4 HSE CG CPH1 6 1.8000 0.0500
12.0110 0.2200 1.7000 0.7200
58 4 HSE CE1 CPH2 6 1.8000 0.0500
12.0110 0.2500 1.7000 0.7200
59 4 HSE HE1 HR1 1 0.9000 0.0460
1.0080 0.1300 1.3000 0.8500
60 4 HSE NE2 NR1 1 1.8500 0.2000
14.0070 -0.3600 1.3000 0.8500
61 4 HSE HE2 H 1 0.2245 0.0460
1.0080 0.3200 1.2000 0.8500
62 4 HSE CD2 CPH1 6 1.8000 0.0500
12.0110 -0.0500 1.7000 0.7200
63 4 HSE HD2 HR3 1 1.4680 0.0078
1.0080 0.0900 1.3000 0.8500
64 4 HSE C C 6 2.0000 0.1100
12.0110 0.5100 1.7000 0.7200
65 4 HSE O O 8 1.7000 0.1200
15.9990 -0.5100 1.5000 0.8500
Note atom #56 -- ND1. The atomic number is listed as 1 (from the
ATOMIC_NUMBER field). If this comes from chamber in ParmEd, then there's
probably something wrong with the RTF file defining the HSE residue (is
that a Histidine protomer?). The atomic number is taken from the MASS
section of either the parameter or RTF files, falling back to guessing it
from the atomic mass if the element name is not specified. Since the mass
is that of Nitrogen, the element name must be set as H in those files :(.
I'd need to see the original files where this prmtop was generated from in
order to find out why ND1 and NR1 are being treated like hydrogens. But
this would certainly explain what you were seeing.
In the prmtop you attached, I found the following atoms that ParmEd
identified as "hydrogen" (from the ATOMIC_NUMBER flag) that had more than 1
bond partner:
[<Atom ND1 [55]; In HSE 3>, <Atom NE2 [59]; In HSE 3>, <Atom ND1 [535]; In
HSE 36>, <Atom NE2 [539]; In HSE 36>, <Atom ND1 [891]; In HSE 58>, <Atom
NE2 [895]; In HSE 58>, <Atom ND1 [1729]; In HSE 113>, <Atom NE2 [1733]; In
HSE 113>, <Atom ND1 [1812]; In HSE 118>, <Atom NE2 [1816]; In HSE 118>,
<Atom ND1 [2219]; In HSE 147>, <Atom NE2 [2223]; In HSE 147>, <Atom ND1
[2386]; In HSE 159>, <Atom NE2 [2390]; In HSE 159>, <Atom ND1 [2993]; In
HSE 197>, <Atom NE2 [2997]; In HSE 197>, <Atom ND1 [4369]; In HSE 291>,
<Atom NE2 [4373]; In HSE 291>, <Atom ND1 [4386]; In HSE 292>, <Atom NE2
[4390]; In HSE 292>, <Atom ND1 [4841]; In HSE 320>, <Atom NE2 [4845]; In
HSE 320>, <Atom ND1 [5272]; In HSD 350>, <Atom NE2 [5277]; In HSD 350>,
<Atom ND1 [5587]; In HSE 368>, <Atom NE2 [5591]; In HSE 368>, <Atom ND1
[5916]; In HSE 388>, <Atom NE2 [5920]; In HSE 388>, <Atom ND1 [6214]; In
HSE 408>, <Atom NE2 [6218]; In HSE 408>, <Atom ND1 [6231]; In HSE 409>,
<Atom NE2 [6235]; In HSE 409>, <Atom ND1 [6284]; In HSE 412>, <Atom NE2
[6288]; In HSE 412>, <Atom ND1 [6678]; In HSD 438>, <Atom NE2 [6683]; In
HSD 438>]
It looks like the 2 histidine protomers are the culprits here. It also
looks like the shakeH message is being printed correctly by pmemd.cuda, but
there's something wrong in the source CHARMM topology file.
Thanks,
Jason
--
Jason M. Swails
_______________________________________________
AMBER-Developers mailing list
AMBER-Developers.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber-developers
Received on Tue May 22 2018 - 06:30:02 PDT
