********> bugfix.68
Author:   Jim Caldwell
Date:     3/22/96

Programs: Parm

Severity: Extreme

Problem:  When setting up pert parm files, the use of "-m frcmod" 
          to enter any specific dihedral terms may ruin some or
          all of the dihedrals in the "pert" group.

          To discover if this bug has affected a given perturbation,
          compare the dihedral energies for the initial step of Sander 
          and Gibbs using the same input coordinate file. 

Note:     Leap does not have this bug.

Workaround: Use Leap instead of Parm or enter new force field terms in 
          the "-f parm.dat" force field file.

History:  Amber has always had the dangerous restriction that
          dihedral parameters must be entered as 'general' types 
          (X-a-b-X, where 'X' can be any atom) first, followed 
          by the 'specific' types (c-a-b-d).  The 'normal', non-
          perturbation part of parm was fixed by Release4.0/bugfix.46,
          but the area dealing with perturbed groups wasn't.

          The new code sorts the dihedrals so that they may be entered
          in any order via -f and/or -m files. 

          Note: this changes the order of terms in the prmtop file, 
          so amber41/test/parm/ diffs with the demo tree issued before 
          this fix will appear to fail; the results have been tested by 
          copying the new prmtop files back into the demo tree and
          rerunning sander and gibbs tests to verify that the energies
          were the same.

Fix:  Make the following changes to parmr.f
----------------------------------------------------------------------------
*** OLD parmr.f
--- NEW parmr.f
***************
*** 61,66 ****
--- 61,72 ----
  C
        DIMENSION JSOLTY(MAXTYP),ISOLTY(MAXTYP),MBT(4)
        DIMENSION IXCHG(MAXTYP),JXCHG(MAXTYP)
+       dimension is_temp(maxdpt),js_temp(maxdpt),ks_temp(maxdpt),
+      -          ls_temp(maxdpt)
+       dimension pks(3*maxdpt),phases(3*maxdpt),pns(3*maxdpt)
+       dimension ix_temp(maxdpt),jx_temp(maxdpt),kx_temp(maxdpt),
+      -          lx_temp(maxdpt)
+       dimension pkx(3*maxdpt),phasex(3*maxdpt),pnx(3*maxdpt)
  c
  c     Rev A mod:  use double temps
        double precision atmwt,trk,treq,ttk,tteq,tpk,tphase,tpn,ttpn,a,b
***************
*** 482,487 ****
--- 488,599 ----
        READ(in2,18) TPK,TPHASE,TPN
        GO TO 1313
    107 CONTINUE
+ 
+ c now sort the diehdrals.
+       kn = 1
+       ls = 1
+       kx = 1
+       lx = 1
+ 
+       ms = 1
+       mx = 1
+       do jn = 1,i
+          if(ipt(jn) .ne. iatx) then
+            is_temp(ls) = ipt(jn)           
+            js_temp(ls) = jpt(jn) 
+            ks_temp(ls) = kpt(jn) 
+            ls_temp(ls) = lpt(jn) 
+            ls = ls  + 1
+ 777        if (pn(kn)  .lt. 0.0)  then
+               pks(ms) = pk(kn) 
+               phases(ms) = phase(kn)
+               pns(ms) = pn(kn)
+               ms = ms + 1
+               kn = kn + 1
+               go to 777
+            endif
+ c                                   put the one pn + in
+            pks(ms) = pk(kn) 
+            phases(ms) = phase(kn)
+            pns(ms) = pn(kn)
+            ms = ms + 1
+            kn = kn + 1
+         else
+            ix_temp(lx) = ipt(jn)           
+            jx_temp(lx) = jpt(jn) 
+            kx_temp(lx) = kpt(jn) 
+            lx_temp(lx) = lpt(jn) 
+            lx = lx  + 1
+ 7777       if (pn(kn)  .lt. 0.0)  then
+               pkx(mx) = pk(kn) 
+               phasex(mx) = phase(kn)
+               pnx(mx) = pn(kn)
+               mx = mx + 1
+               kn = kn + 1
+               go to 7777
+            endif
+ c                                   put the one pn + in
+            pkx(mx) = pk(kn) 
+            phasex(mx) = phase(kn)
+            pnx(mx) = pn(kn)
+            mx = mx + 1
+            kn = kn + 1
+         endif
+       enddo
+ c                                     now pack the dihedrals back
+       lim = lx-1
+       kn = 1
+       mn = 1
+       ln = 1
+       do jn = 1,lim
+            ipt(kn) = ix_temp(jn)
+            jpt(kn) = jx_temp(jn)
+            kpt(kn) = kx_temp(jn)
+            lpt(kn) = lx_temp(jn)
+            ipntr(kn) = mn 
+            kn = kn + 1
+ 888        if(pnx(ln)  .lt. 0.0)then
+               pk(mn) = pkx(ln) 
+               phase(mn) = phasex(ln)
+               pn(mn) = pnx(ln)
+               mn = mn + 1
+               ln = ln + 1
+               go to 888
+            endif
+            pk(mn) = pkx(ln) 
+            phase(mn) = phasex(ln)
+            pn(mn) = pnx(ln)
+            mn = mn + 1
+            ln = ln + 1
+       enddo
+ 
+       lim = ls-1
+ 
+       ln = 1
+       do jn = 1,lim
+               ipt(kn) = is_temp(jn)
+               jpt(kn) = js_temp(jn)
+               kpt(kn) = ks_temp(jn)
+               lpt(kn) = ls_temp(jn)
+               ipntr(kn) = mn 
+               kn = kn + 1
+ 1888          if(pns(ln)  .lt. 0.0)then
+                  pk(mn) = pks(ln) 
+                  phase(mn) = phases(ln)
+                  pn(mn) = pns(ln)
+                  mn = mn + 1
+                  ln = ln + 1
+                  go to 1888
+               endif
+               pk(mn) = pks(ln) 
+               phase(mn) = phases(ln)
+               pn(mn) = pns(ln)
+               mn = mn + 1
+               ln = ln + 1
+       enddo
+ 
+       ipntr(i+1)=mn
+ 
        if (iparml .ge. 1) write(iout,'(/12x,a,/)') 'Improper Dihedrals:'
        NUMPHI = I
        NPTRA = IPNTR(I+1)-1
----------------------------------------------------------------------------

--