AMBER Archive (2003)

Subject: AMBER: mm_pbsa - ?small bugs?

From: Chris Moth (Chris.Moth_at_vanderbilt.edu)
Date: Wed Sep 24 2003 - 12:37:37 CDT

I am taking my first steps with MM_PBSA and I think I have encountered some
small bugs in the script and make_crd_hg.f files. Apologies in advance if
these have already been reported or are my errors!

At our site, the executable to build trajectories is called make_crd_hg -
not make_crd. Until I changed

@GC_FIL = ("MAKE_CRD");

to

@GC_FIL = ("MAKE_CRD_HG");

in mm_pbsa.pl, I was unable to generate trajectory snapshots due to the
message:

        "File MAKE_CRD is not found"

I made this substitution three other places in the .pl script. Please let
me know if this is a mistake on my part, or if we have built the mm_pbsa
suite incorrectly here.

================================================================
I believe there may also be a bug in make_crd_hg.f

I'll use the 01_Example to demonstrate this - and avoid getting into
details of my own work.

If I change mm_pbsa.in from

NUMBER_REC_GROUPS 2
RSTART 1
RSTOP 2621
RSTART 3863
RSTOP 3907

to

NUMBER_REC_GROUPS 1
RSTART 1
RSTOP 2621

I get lots of happy messages followed by:

=>> Creating input
     make_crd input

=>> Creating coordinates
     Executing makecrd
  Something wrong with NGR 2622 2621
Broken pipe

I think my mm_pbsa.in change should be "legal" If not, please let me know
why it is not. Otherwise, read on.

I have checked the make_crd_hg.f source code and I believe the error may be
in these statements:

           j=1
           k=1
           NGR=0
           NGL=0
           NGC=0
           do i=1,NGREC+NGLIG
             if(((NGREC.gt.0).and.(RSTART(j).lt.LSTART(k))).or.
      1 (k.gt.NGLIG)) then
               do l=RSTART(j),RSTOP(j)
                 if(x(l).lt.big) then
                   NGR=NGR+1
                   xR(NGR)=x(l)
                   yR(NGR)=y(l)
                   zR(NGR)=z(l)
                   if(NGCOM.gt.0) then
                     NGC=NGC+1
                     xC(NGC)=x(l)
                     yC(NGC)=y(l)
                     zC(NGC)=z(l)
                   endif
                 endif
               end do
               j=j+1
             elseif(((NGLIG.gt.0).and.(LSTART(k).lt.RSTART(j))).or.
      1 (j.gt.NGREC)) then
               do l=LSTART(k),LSTOP(k)
                 if(x(l).lt.big) then
                   NGL=NGL+1
                   xL(NGL)=x(l)
                   yL(NGL)=y(l)
                   zL(NGL)=z(l)
                   if(NGCOM.gt.0) then
                     NGC=NGC+1
                     xC(NGC)=x(l)
                     yC(NGC)=y(l)
                     zC(NGC)=z(l)
                   endif
                 endif
               end do
               k=k+1
             endif
           end do

C
C Check again
C

========================================================
I don't understand all the nuances of FORTRAN :) - but I think structurally
the error could be avoided with two separate
loops over NGR and NGL so that the dependence of the variables on both
NGREC and NGLIG (number of receptor and ligand groups)
could be separated out. Also, I think we can get rid of j and k:

Something like this might do the trick - though I have not compiled/tested
this:

           NGR=0
           NGL=0
           NGC=0
           if (NGREC.gt.0)
           do i=1,NGREC
               do l=RSTART(i),RSTOP(i)
                 if(x(l).lt.big) then <<<<< Let me add - if x(l).nlt.big
don't we have an error? - and shouldn't the program exit with a message?
                   NGR=NGR+1
                   xR(NGR)=x(l)
                   yR(NGR)=y(l)
                   zR(NGR)=z(l)
                   if(NGCOM.gt.0) then
                     NGC=NGC+1
                     xC(NGC)=x(l)
                     yC(NGC)=y(l)
                     zC(NGC)=z(l)
                   endif
               end do
            end do

           if (NGLIG.gt.0)
           do i=1,NGLIG
               do l=LSTART(i),LSTOP(i)
                 if(x(l).lt.big) then <<<<< Let me add - if x(l).nlt.big
don't we have an error? - and shouldn't the program exit with a message?
                   NGL=NGL+1
                   xL(NGL)=x(l)
                   yL(NGL)=y(l)
                   zL(NGL)=z(l)
                   if(NGCOM.gt.0) then
                     NGC=NGC+1
                     xC(NGC)=x(l)
                     yC(NGC)=y(l)
                     zC(NGC)=z(l)
                   endif
                 endif
               end do
           end do
        endif

C
C Check again
C

Thanks

Chris Moth
http://www.structbio.vanderbilt.edu/~cmoth/mddisplay

-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber_at_scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo_at_scripps.edu