> On Oct 31, 2014, at 4:36 PM, Pete Kekenes-Huskey <huskeypm.gmail.com> wrote:
>
> Hi Jason,
> Thanks for your continued help. We have made progress with your
> suggestions, but are still confronting problems with SHAKE. Namely,
> minimization/equilibration appear to run fine with pmemd (mpi), but the
> same equilibration procedure fails with pmemd.cuda (Amber 14):
>
> 785: 31632 - 31633 and 787: 31633 - 31634
>
> Hydrogen atom 31633 appears to have multiple bonds to atoms 31632 and 31634
> which is illegal for SHAKEH.
>
> Exiting due to the presence of inconsistent SHAKEH hydrogen clusters.
>
> I've verified that the guilty hydrogen should indeed be bonded to both the
> oxygen and opposing hydrogen (vmd output below) and that the geometry
> appears normal. Oddly, this error is thrown for the TIP3 water with resid
> 10001, so there are 10K TIP3 waters before that entry that appear to be ok.
> The .in file we ran was the following:
>
> &cntrl
>
> imin = 0,
>
> irest = 0,
>
> ntx = 1,
>
> ntb = 1,
>
> cut = 10,
>
> ntr = 1,
>
> ntc = 2,
>
> ntf = 2,
>
> ig=-1,
>
> tempi = 0.0,
>
> temp0 = 300.0,
>
> ntt = 3,
>
> gamma_ln = 1.0,
>
> nstlim = 1000, dt = 0.002,
>
> watnam='TIP3',OWTNM='OH2',
>
> hwtnm1='H1',hwtnm2='H2',
>
> ntpr = 100, ntwx = 100, ntwr = 1000
>
> /
>
> Hold the protein fixed with weak restraints
>
> 10.0
>
> RES 1 90
>
> END
>
> END
>
> /////
>
> I've uploaded all the files we've used to
> http://users.mccammon.ucsd.edu/~huskeypm/share/filez.tgz. To summarize, we
> created the input files with CHARMM-GUI (in order to be consistent with
> previous runs) and the charmm-36 FF. We used parmed (latest version from
> the git repo) to obtain prmtop/inpcrd files:
>
> chamber -top toppar/top_all36_prot.rtf -param toppar/par_all36_prot.prm
> -str toppar/toppar_water_ions.str -psf step2_solvator.xplor.psf -crd
> step2_solvator.crd -box 80.,80.,80.
>
> outparm apo.prmtop apo.inpcrd
>
> Minimizations were done using pmemd. All input files should be in the link
> I provided above.
>
> If you're able to look into this and get this running via pmemd.cuda, we'd
> be very grateful. I'm sure were missing something subtle, so a third pair
> of eyes would be helpful.
OK, so this was a very weird problem. The problem is with step2_solvator.xplor.psf. This PSF file is not “extended” (note the lack of EXT in the top line, which is present in the other step2_*.psf files), which means the residue number field can have at most 4 characters (so no chain can have more than 9999 residues). Look around line 31637 of step2_solvator.xplor.psf, where it appears that the last digit on the residue number is simply shaved off (so that 10 water residues are conflated into the single residue 10000).
As a result, ParmEd puts these 30 atoms into a single residue (since that’s what the PSF files tells it to do). To pmemd.cuda, this does not look like a water (that should be single residues of exactly 3 atoms). As a result, those constraints do not get dispatched to the fast SETTLE routines and the “slow” (i.e., approximate, iterative SHAKE) routine in pmemd.cuda is _not_ generalized to handle hydrogen bonded to more than one atom.
The solution is to use literally _either_ of the other two step2_solvator PSF files (which are extended). :) I tested it with both options and it solved the problem you were seeing on my machine with pmemd.cuda.
HTH,
Jason
--
Jason M. Swails
BioMaPS,
Rutgers University
Postdoctoral Researcher
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Sat Nov 01 2014 - 16:00:02 PDT
