Dave,
one thing to keep in mind is the limitations on shake if
you're using MPI. do the constrained atoms have any H attached
to them? Are they in different residues? That won't work due to
the iterative way shake is done. Try single thread if you haven't
already. Also you could put some debugging into the bond force
calculation to see if it actually uses your force constant.
try using rdparm to print out the bond lists in your new and old prmtop.
carlos
On Fri, Feb 22, 2008 at 10:15 PM, David Cerutti
<dcerutti.mccammon.ucsd.edu> wrote:
> So, I've now got a program that seems to be able to parse an AMBER prmtop
> file and print out a new one, with at least one new bond, that will be
> properly read and used by AMBER. A new bond is defined by:
>
> 1.) Dealing with the bonds:
> - Incrementing NUMBND by 1.
> - Incrementing NBONH by 1.
> - Adding the stiffness and equilibrium length of the new bond
> to the BOND_FORCE_CONSTANT and BOND_EQUIL_VALUE fields, respectively. I
> just use 100 kcal/mol-A^2 for the stiffness, because I don't expect it to
> matter if the bond length is to be constrained. The new values for
> stiffness and length are added at the end of their respective parameter
> lists, so the new bond index number is the new value of NUMBND.
> - Adding the new bond to the very end of the BONDS_INC_HYDROGEN list.
> Inserting a bond between atoms 3 and 17 for a topology that ORIGINALLY had
> 39 bond stiffness/equilibrium values (and now has 40 including my new
> bond) would imply I add " 9 51 40" to the end of the
> BONDS_INC_HYDROGEN list.
>
> 2.) Dealing with the exclusions:
> - Making space for the new exclusion in the excluded atoms list and
> catalogging it there.
> - Incrementing the number of exclusions for the lower-numbered atom in the
> bond (i.e. adding a bond between atoms 3 and 17 implies I increment the
> number of exclusions for atom 3 by 1).
> - Incrementing the total number of exclusions by 1.
>
> As I said, the resulting file can be read by AMBER and it appears that a
> new bond I introduce in this way have some effect. However, the new bond
> is not being SHAKE-constrained like I want it to be. Since I have added
> this bond to the "BONDS_INC_HYDROGEN" field I am running with ntf = 2, ntc
> = 2 to turn on SHAKE for bonds containing hydrogen. Apparently, however,
> this is not sufficient. Another thing that I'm confused about is that the
> stiffness I enter for the new bond does not appear to have any effect, but
> the equilibrium length I enter does--if I insert a bond between two atoms
> and then run some MD and monitor the distance between those atoms, I can
> compare it to the result obtained with the original topology. Doing so, I
> find I can pull those atoms closer together or push them further apart by
> inserting a bond that is shorter or longer than the initial distance
> between the atoms.
>
> I think I'm really close to having this working. Hoping anyone can point
> out where I may have gone wrong...
>
> Thanks,
>
> Dave
> -----------------------------------------------------------------------
> The AMBER Mail Reflector
> To post, send mail to amber.scripps.edu
> To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
>
-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber.scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
Received on Sun Feb 24 2008 - 06:07:54 PST
