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
