Dear Dr. Case,
I wanted to thank you again for your advice on adding new atom types. I took a quick peek at the source of program Leap and ParmED, and I think the PB radii assignment seems dependent on the atomic number. Thus, it’s necessary to make sure that the atomic number is correct, even though in my case it doesn’t do anything to the other parameters for the simulation.
I have also figured out the previous problem I had with my prmtop file. The cause is likely a bug in parmchk2 for not commenting all ATTN parameters in the dihedral section, but only the first occurrence.
Please see below as part of my frcmod file:
CT-CT-OH-HO 1 0.160 0.000 -3.000
CT-CT-OH-HO 1 0.250 0.000 1.000
CT-CT-SO-O 1 0.000 0.000 2.000 ATTN, need revision
CT-CT-SO-O 0 0.000 0.000 0.000
H1-CT-CT-OH 1 0.000 0.000 -3.000
H1-CT-CT-OH 1 0.250 0.000 1.000
H1-CT-CT-SO 9 1.400 0.000 3.000
H1-CT-CT-H2 9 1.400 0.000 3.000
CT-CT-CT-OH 9 1.400 0.000 3.000
CT-CT-CT-SO 9 1.400 0.000 3.000
CT-CT-CT-H2 9 1.400 0.000 3.000
OH-CT-SO-O 1 0.000 0.000 2.000 ATTN, need revision
OH-CT-SO-O 0 0.000 0.000 0.000
SO-CT-OH-HO 3 0.500 0.000 3.000
N -C -CT-HC 6 0.000 0.000 2.000
O -C -CT-HC 1 0.800 0.000 -1.000
O -C -CT-HC 1 0.000 0.000 -2.000
O -C -CT-HC 1 0.080 180.000 3.000
H2-CT-OH-HO 3 0.500 0.000 3.000
H2-CT-SO-O 1 0.000 0.000 2.000 ATTN, need revision
H2-CT-SO-O 0 0.000 0.000 0.000
The uncommented dihedral parameters with all zeros could then leak to the prmtop file… The solution is very simple in my case: Remove all parameters about SO from the parm10 frcmod. I should have checked the frcmod files more carefully, before using the files for high throughput calculations…
I’m using Amber 20. Not sure if this minor issue has been fixed in the recent version. I can open a separate topic and provide more details on this if you think that could be of some help. Thank you again for your time and kind advice!
Bests,
Amy
From: He, Amy via AMBER <amber.ambermd.org>
Date: Monday, September 5, 2022 at 1:21 PM
To: David A Case <david.case.rutgers.edu>, AMBER Mailing List <amber.ambermd.org>
Subject: Re: [AMBER] ATOMIC NUMBER section in prmtop file - sulfur atom in ligand is -1?
Dear Dr. Case,
Thank you so much for your thoughtful reply! That was really helpful. I thought the problem was relevant to hydrogen mass partitioning, but as you correctly pointed out, it is likely more about the parametrization of the sulfur atom.
But I’m still not so sure what to do about the problem… And I would really appreciate it if you could give some more advice on it. My apologies for the lengthy message in advance.
1. New atom type…?
I did notice that the ligand bisulfite sulfur could be a new atom type. The sulfur has atom type SO in my prepin file, as assigned from Antchamber, which is not in $AMBERHOME/dat/leap/parm/parm10.dat, but in $AMBERHOME/dat/antechamber/ATOMTYPE_AMBER.DEF.
Therefore, Antechamber didn’t reject the atom type. In addition, I was able to use parmchk2 so that the relevant parameters are complemented with the GAFF2 equivalents.
However, without declaring the new atom type in Leap, there would be a message in the standard output:
# +Currently only Sp3-Sp3/Sp3-Sp2/Sp2-Sp2 are supported
# +---Tried to superimpose torsions for: *-C55-S58-*
# +--- With Sp3 - Sp0
# +--- Sp0 probably means a new atom type is involved
# +--- which needs to be added via addAtomTypes
Still, Leap was able to write a prmtop file for my system, though it seems like it could not recognize the hybridization state of the sulfur… I disregarded this message the first time I ran the calculation, because I thought this torsion had been taken care of in the GAFF2.frcmod I made:
H2-CT-SO-O 9 1.300 0.000 3.000 same as X -c3-s6-X , penalty score= 0.0
2. My Leap input file
cat << EOF > leap.addatomtype.in
source leaprc.protein.ff14SB
source leaprc.gaff2
source leaprc.water.tip3p
addAtomTypes {
{ "SO" "S" "sp3" }
}
set default PBRadii mbondi3
loadamberprep RES.prepin
loadamberparams RES_gaff2.frcmod
loadamberparams RES_ff14SB.frcmod
complex = loadpdb complex.pdb4amber.pdb
saveamberparm complex prmtop.prmtop inpcrd.inpcrd
quit
EOF
3. Difference in prmtop files
I used a compare tool to examine the prmtop files. I found that addAtomTypes changed ATOMIC_NUMBER, RADII, and SCREEN in the prmtop.
ATOMIC_NUMBER: from -1 to 16
RADII: from 1.5 to 1.8
SCREEN: from 0.8 to 0.96
4. Results of the minimization
I found a big difference in the EEL term, with or without adding the new atom type in Leap!!!
(1) No Leap addAtomTypes, No ParmED HMR
NSTEP ENERGY RMS GMAX NAME NUMBER
1000 -1.0844E+04 9.4968E-02 1.5259E+00 CE2 3311
BOND = 164.5286 ANGLE = 594.5273 DIHED = 3359.7041
VDWAALS = -2761.2412 EEL = -21251.9696 EGB = -4408.7354
1-4 VDW = 934.7608 1-4 EEL = 11987.3438 RESTRAINT = 0.0000
ESURF = 536.9040
(2) Leap addAtomTypes, No ParmED HMR
NSTEP ENERGY RMS GMAX NAME NUMBER
1000 -7.2047E+04 2.9975E+00 1.0349E+02 O60 4654
BOND = 316.5836 ANGLE = 1018.8918 DIHED = 3472.6794
VDWAALS = -2755.0830 EEL = -83097.9364 EGB = -4361.2372
1-4 VDW = 865.3401 1-4 EEL = 11956.3905 RESTRAINT = 0.0000
ESURF = 537.265
With adding the new atom type in Leap, the endpoint structure appeared to be really distorted, just like it was with the ParmED prmtop file. And I thought that was not just because of the changes in RADII and SCREEN.
So I was kind of thinking, maybe ATOMIC NUMBER is -1 is the correct option (and it is necessary) for the sulfur atom, because it doesn’t have a type in parm10, and has been overridden by its equivalent in GAFF2?? Or maybe I did something that was serious wrong in the parameterization :’(
I’m sharing my input files and commands:
Ligand pdb as input for Antechamber
https://urldefense.com/v3/__https://drive.google.com/file/d/1cTxnJdA9a4NHLH4NRZfUNRXdvu-T6p3L/view?usp=sharing__;!!KGKeukY!3gF2ENRQYw2tQS0MTpggbBUm2W6_UbSz9zF_gMA3leVPDs9LwWZd-IRal9COzll9OVyb7RKt-n8AIWOdHuKT2eKJ-_Y$<
https://urldefense.com/v3/__https:/drive.google.com/file/d/1cTxnJdA9a4NHLH4NRZfUNRXdvu-T6p3L/view?usp=sharing__;!!KGKeukY!3gF2ENRQYw2tQS0MTpggbBUm2W6_UbSz9zF_gMA3leVPDs9LwWZd-IRal9COzll9OVyb7RKt-n8AIWOdHuKT2eKJ-_Y$>
Complex pdb as input for Leap
https://urldefense.com/v3/__https://drive.google.com/file/d/1-Ya_aL4vhvvg--umDozso_GUhoiNli6Y/view?usp=sharing__;!!KGKeukY!3gF2ENRQYw2tQS0MTpggbBUm2W6_UbSz9zF_gMA3leVPDs9LwWZd-IRal9COzll9OVyb7RKt-n8AIWOdHuKTeXNVrFQ$<
https://urldefense.com/v3/__https:/drive.google.com/file/d/1-Ya_aL4vhvvg--umDozso_GUhoiNli6Y/view?usp=sharing__;!!KGKeukY!3gF2ENRQYw2tQS0MTpggbBUm2W6_UbSz9zF_gMA3leVPDs9LwWZd-IRal9COzll9OVyb7RKt-n8AIWOdHuKTeXNVrFQ$>
Bash commands I used to execute programs and generate other inputs
https://urldefense.com/v3/__https://drive.google.com/file/d/1cTxnJdA9a4NHLH4NRZfUNRXdvu-T6p3L/view?usp=sharing__;!!KGKeukY!3gF2ENRQYw2tQS0MTpggbBUm2W6_UbSz9zF_gMA3leVPDs9LwWZd-IRal9COzll9OVyb7RKt-n8AIWOdHuKT2eKJ-_Y$<
https://urldefense.com/v3/__https:/drive.google.com/file/d/1cTxnJdA9a4NHLH4NRZfUNRXdvu-T6p3L/view?usp=sharing__;!!KGKeukY!3gF2ENRQYw2tQS0MTpggbBUm2W6_UbSz9zF_gMA3leVPDs9LwWZd-IRal9COzll9OVyb7RKt-n8AIWOdHuKT2eKJ-_Y$>
I hope the formatting of this message is readable.. The outcome should be reproducible with the shared files. Once again, massive thanks to you for your time and kind advice.
Many Thanks,
Amy
--
Amy He
Chemistry Graduate Teaching Assistant
Hadad Research Group
Ohio State University
he.1768.osu.edu
From: David A Case <david.case.rutgers.edu>
Date: Monday, September 5, 2022 at 9:14 AM
To: He, Amy <he.1768.buckeyemail.osu.edu>, AMBER Mailing List <amber.ambermd.org>
Subject: Re: [AMBER] ATOMIC NUMBER section in prmtop file - sulfur atom in ligand is -1?
On Mon, Sep 05, 2022, He, Amy via AMBER wrote:
>One is the original prmtop I generated by Leap. In the ATOMIC NUMBER
>section, the ligand sulfur atom is -1.
I'm guessing(?) that you defined a new sulfur atom type for the ligand,
but neglected to have a addAtomTypes command in your tleap script:
addAtomTypes {
{ "xx" "S" "sp3" }
}
Where "xx" is your new atom type. Still, this should not make any
difference, since (I think!) Amber only really uses the atomic number for
QM/MM calculations, not for conventional force fields.
Also, the masses are irrelevant for minimization. I think we would need
to have files to reproduce the problem you are seeing (for a small test
case). It sounds like my assumption that the atomic number section is not
relevant is probably wrong.
You might be able to make some progress by looking the breakdown of energy
terms for the starting structure, using both prmtop files. Are all fields
the same except for the VDW term?
If you did leave out the addAtomTypes entry, you could try to fix that, then
regenerate the prmtop file with tleap. Maybe something else is different as
well. Use the parmed commands printLJTypes and printLJMatrix to examine
your nonbonded interactions. Look especially for things related to the
sulfur atom.
...good luck...dac
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
https://urldefense.com/v3/__http://lists.ambermd.org/mailman/listinfo/amber__;!!KGKeukY!3gF2ENRQYw2tQS0MTpggbBUm2W6_UbSz9zF_gMA3leVPDs9LwWZd-IRal9COzll9OVyb7RKt-n8AIWOdHuKT8jrmTzc$<https://urldefense.com/v3/__http:/lists.ambermd.org/mailman/listinfo/amber__;!!KGKeukY!3gF2ENRQYw2tQS0MTpggbBUm2W6_UbSz9zF_gMA3leVPDs9LwWZd-IRal9COzll9OVyb7RKt-n8AIWOdHuKT8jrmTzc$>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Tue Sep 06 2022 - 13:30:02 PDT
