Thank you for your reply.
Yes, I am aware of the fact that VdW parameters are 0. I guess you could call it a prmtop-based FEP scheme, that's not a bad approximation.
As for the hydrogens being in the "wrong" list, this is deliberate. We were having problems with SHAKE interfering with our custom atoms, so we needed those bonds to not be considered as including hydrogens. This was accomplished by setting the element of the custom hydrogens to gadolinium (hence the atomic number of 64), though their mass and other parameters are then set to those of a normal hydrogen. The result is that the atoms are effectively hydrogens but are not thought of as such when writing the prmtop. Perhaps bad practice, but as a workaround it works well enough for the moment.
Is putting them in the "DIHEDRALS INCLUDING HYDROGEN" list the only solution?
Pietro
-----Original Message-----
From: amber-bounces.ambermd.org [mailto:amber-bounces.ambermd.org] On Behalf Of Jason Swails
Sent: 23 July 2013 15:30
To: AMBER Mailing List
Subject: Re: [AMBER] Parmed.py deletedihedral
On Tue, Jul 23, 2013 at 7:43 AM, Aronica, Pietro < pietro.aronica07.imperial.ac.uk> wrote:
> I have sent the topology file.
>
OK, I've found the problem. The issue is in the VXI (custom) residue, as I'm sure you're aware. The problem is that you used the atom names "mh"
and "dh" to denote hydrogen atoms, but you did not tell tleap that those atoms were supposed to be hydrogens. (Indeed, if you look at the ATOMIC_NUMBER section, all of those atom types have atomic number 64.)
As a result, tleap put the dihedral you're looking for (.33 -- .32 -- .36
-- .37) into the DIHEDRALS_WITHOUT_HYDROGEN section of the topology file.
ParmEd, on the other hand, looks at the masses of the atoms, determines that dh and mh types are actually hydrogens, and proceeds to look for them in the DIHEDRALS_INC_HYDROGEN section of the prmtop... which it cannot find.
You should go back to your original leap file that you used to generate the prmtop and add the following command to the beginning of the file:
addAtomTypes {
{ "dh" "H" "sp3" }
{ "mh" "H" "sp3" }
...
}
for all of the atom types that you added to VXI (ones that are naturally in gaff should already be properly loaded inside leaprc.gaff -- you can see an example of the addAtomTypes syntax in there). This is important, since sander and pmemd should know which atoms are Hs and which are not -- if you currently use SHAKE on H-atoms only (ntc == ntf == 2), then sander and pmemd will _not_ apply SHAKE to these atom types as you would expect (at least to my knowledge they won't).
Another issue I found looking at your topology file is that it appears that no LJ parameters were set for several of the atom types
Reading input from STDIN...
> printDetails .%d=,m=,f=
The mask .%d=,m=,f= matches 6 atoms:
ATOM RES RESNAME NAME TYPE LJ Radius LJ Depth Mass
Charge GB Radius GB Screen
32 2 VXI CG2 fc 1.9080 0.1094 12.0100
-0.3192 1.7000 0.7200
33 2 VXI HG21 mh 1.4870 0.0157 1.0080
0.0791 1.2000 0.8500
36 2 VXI CD1 dc 0.0000 0.0000 12.0100
0.0000 1.7000 0.7200
37 2 VXI HD11 dh 0.0000 0.0000 1.0080
0.0000 1.2000 0.8500
38 2 VXI HD12 dh 0.0000 0.0000 1.0080
0.0000 1.2000 0.8500
39 2 VXI HD13 dh 0.0000 0.0000 1.0080
0.0000 1.2000 0.8500
>
If you notice here, the LJ radius and depth for CD1, HD11, HD12, and HD13 are 0. Unless you are specifically testing something here (or disappearing some atoms through a prmtop-based FEP scheme, no C atom should have 0 VDW.
Likewise for H atoms attached to carbons (those attached to oxygens have 0 LJ parameters).
HTH,
Jason
--
Jason M. Swails
BioMaPS,
Rutgers University
Postdoctoral Researcher
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Tue Jul 23 2013 - 08:00:04 PDT
