Re: [AMBER] CHAMBER: error in finding bond type

From: Mark Williamson <mjw.mjw.name>
Date: Wed, 17 Oct 2012 13:29:17 +0100

On 16/10/12 13:24, Mark Williamson wrote:
> On 16/10/12 12:28, Mike Limb wrote:
>> Dear Amber users,
>>
>> I am currently trying to use CHAMBER in order to convert my CHARMM
>> compatible files to use in AMBER.
>>

Dear Mike Limn,

OK, I have a fix and a verbose explanation as to what went wrong; it's
essentially an oversight with a hard limit parameter within the code.



Your input fails with the Error message:

    --- Assigned NONBONDED parameters
  <bond_types>error in finding bond type 1990 0
    0 537
  Cannot continue

 
 

Overview
========

 

This being generated at psfprm.F90:1969 in the bond_types subroutine,
which organises bond types found in the psf file.

.........
    do n=1,nbond
       i=ib(n)
       j=jb(n)
       it=locattype(i)
       jt=locattype(j)
       if(it > jt)call swap(i,j,it,jt)
       icb(n)=ibndt(it,jt)
       if(icb(n) == 0 ) then
          write(outu,*)"<bond_types>","error in finding bond
type",n,it,jt,ioff
          call err_exit("Cannot continue")
       endif
    enddo
............


- The nbond value (22848) is the number of bonds and is set in the
read_psf_bond_angle_dihe routine, called from main:

    if(verbose_debug)write(outu,'(a)')"<read_bond_angle_dihe> bonds"
           read(u,fmt00) nbond

This is it reading from the psf file:

.....
    22848 !NBOND: bonds
        2 1 3 1 4 1 7 5
       10 7 13 10 16 13 19 16
......


- ib() and jb() are lists of i and j'th atom as a function of which bond
they are in i.e. if bond 5 contained atoms 1 and 2, ib(5) would be 1 and
jb(5) would be 2.


- locattype !Map between natom value ---> attypes()
                 ! \--> attype_name()
                 ! \-> element()
                 !Runs over natom

- icb(nbond) is filled by ibndtyp(it,jt) i.e.bond type for types it,jt


Cause
=====
icb(1990) is zero. This was set by ibndt(it,jt), where ibndt(0,0) = 0
So, why are it and jt 0?

 From hewl_mgm1_solv_charmm.psf, i,j are 2000, 2002, which are these atoms:

     2000 MGM1 1 NAM C1 403 0.200000 12.0110 0
   0.00000 -0.301140E-02

     2002 MGM1 1 NAM O1 404 -0.660000 15.9994 0
   0.00000 -0.301140E-02

Why is locattype(2000) and locattype(2002) returning zero?
Looking at locattype via <atom_types> in the (-verbose) output,
it seems there are no entries for atoms 2000 and 2002:

<atom_types> atom 1997 chmtype 400 local type 40
<atom_types> atom 2001 chmtype 400 local type 40
<atom_types> atom 2005 chmtype 400 local type 40
<atom_types> atom 2008 chmtype 400 local type 40

So, it looks as if something has gone wrong during the filling of the
locattype() array in the get_atom_parameters subroutine. The section in
psfprm.F90 ~ line 328:

         if(itype == attypes(n))then
           nfound=nfound+1
           attype_name(n)(1:idlen)=word_type(1:idLen)
           ........

is never entered because the attypes array does not contain the charm
atom type i.e. the integer (401) in nag_nam_protein_top_CMAP.inp:

.....
MASS 400 HAS 1.00800 H ! sugar aliphatic hydrogen TYPE NO. CHANGED
MASS 401 HOS 1.00800 H ! sugar hydroxyl hydrogen TYPE NO. CHANGED
.....

This in turn is because attypes() is never populated with the value in
the atom_types subroutine due to the hardlimit of maxattypes=400
on the att() array within that subroutine.

Solution
========
In your case, increasing maxattypes constant and recompiling should fix
this.




Also, due to the lack of write access, here's a patch to catch this
situation in future.


diff --git a/AmberTools/src/chamber/psfprm.F90
b/AmberTools/src/chamber/psfprm.F90
index c9d9062..92e9ba8 100644
--- a/AmberTools/src/chamber/psfprm.F90
+++ b/AmberTools/src/chamber/psfprm.F90
.@ -293,6 +293,15 @@ module psfprm
             !OR
             !MASS 1 H 1.00800 H
             call next_int(line,linlen,itype)
+
+ ! Check here to see if itype > maxattypes
+ if ( itype > maxattypes) then
+ write(outu,'(a,i5,a,i5)'),"chmtype", itype, &
+ " is greater than maxattypes hardlimit of ",maxattypes
+ write(outu,'(a)'),"Recompile with a larger value if
maxattypes"
+ call err_exit("Cannot continue")
+ endif
+
             call next_word(line,wlen,linlen,word_type)
             if(word_type(1:idlen) ==
tip3_htype_word(1:idlen))tip3_htype=itype
             if(word_type(1:idlen) ==
tip3_otype_word(1:idlen))tip3_otype=itype



P.S. there is still something quirky with the resulting output; i.e. the
molecules NAG and NAM are corrupted and the CYS residues are
crosslinked, yet still have hydrogens attached. This I think is due to
something being wrong with your inputs... I can't get them to display
correctly in VMD. Can you check that they are correct?


Regards,

Mark

_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Wed Oct 17 2012 - 05:30:04 PDT
Custom Search