Fixed energy differs in WHAM-M (EBE).
[unres.git] / source / wham / src-M / energy_p_new.F
index f9dbc12..cac8518 100644 (file)
@@ -1916,10 +1916,15 @@ C 12/26/95 - for the evaluation of multi-body H-bonding interactions
           ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
           ees=ees+eesij
           evdw1=evdw1+evdwij
-cd          write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
-cd     &      iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
-cd     &      1.0D0/dsqrt(rrmij),evdwij,eesij,
-cd     &      xmedi,ymedi,zmedi,xj,yj,zj
+c             write (iout,'(a6,2i5,0pf7.3,2i5,2e11.3)') 
+c     &'evdw1',i,j,evdwij
+c     &,iteli,itelj,aaa,evdw1
+
+c              write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
+c          write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
+c     &      iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
+c     &      1.0D0/dsqrt(rrmij),evdwij,eesij,
+c     &      xmedi,ymedi,zmedi,xj,yj,zj
 C
 C Calculate contributions to the Cartesian gradient.
 C
@@ -2265,8 +2270,10 @@ C Check the loc-el terms by numerical integration
 C Contribution to the local-electrostatic energy coming from the i-j pair
           eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
      &     +a33*muij(4)
-cd          write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
-cd          write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3)
+c          write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
+c          write (iout,'(a6,2i5,0pf7.3)')
+c     &            'eelloc',i,j,eel_loc_ij
+c          write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3)
           eel_loc=eel_loc+eel_loc_ij
 C Partial derivatives in virtual-bond dihedral angles gamma
           if (calc_grad) then
@@ -2835,7 +2842,9 @@ C Uncomment following three lines for Ca-p interactions
             evdw2_14=evdw2_14+e1+e2
           endif
           evdwij=e1+e2
-c          write (iout,*) i,j,evdwij
+c          write (iout,'(a6,2i5,0pf7.3,2i3,3e11.3)')
+c     &        'evdw2',i,j,evdwij,iteli,itypj,fac,aad(itypj,iteli),
+c     &       bad(itypj,iteli)
           evdw2=evdw2+evdwij
           if (calc_grad) then
 C
@@ -3170,8 +3179,8 @@ C
      & delthe0,sig0inv,sigtc,sigsqtc,delthec,it
       double precision y(2),z(2)
       delta=0.02d0*pi
-      time11=dexp(-2*time)
-      time12=1.0d0
+c      time11=dexp(-2*time)
+c      time12=1.0d0
       etheta=0.0D0
 c      write (iout,*) "nres",nres
 c     write (*,'(a,i2)') 'EBEND ICG=',icg
@@ -3197,8 +3206,8 @@ C Zero the energy function and its derivative at 0 or pi.
         if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
 #ifdef OSF
           phii=phi(i)
-          icrc=0
-          call proc_proc(phii,icrc)
+c          icrc=0
+c          call proc_proc(phii,icrc)
           if (icrc.eq.1) phii=150.0
 #else
           phii=phi(i)
@@ -3212,8 +3221,8 @@ C Zero the energy function and its derivative at 0 or pi.
         if (i.lt.nres .and. itype(i).ne.ntyp1) then
 #ifdef OSF
           phii1=phi(i+1)
-          icrc=0
-          call proc_proc(phii1,icrc)
+c          icrc=0
+c          call proc_proc(phii1,icrc)
           if (icrc.eq.1) phii1=150.0
           phii1=pinorm(phii1)
           z(1)=cos(phii1)
@@ -3281,7 +3290,7 @@ c     &    rad2deg*phii,rad2deg*phii1,ethetai
         if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1
         if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2
         gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)
- 1215   continue
+c 1215   continue
       enddo
 C Ufff.... We've done all this!!! 
       return
@@ -3424,7 +3433,9 @@ C
       etheta=0.0D0
 c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
       do i=ithet_start,ithet_end
-        if (itype(i-1).eq.ntyp1) cycle
+c        if (itype(i-1).eq.ntyp1) cycle
+        if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or.
+     &(itype(i).eq.ntyp1)) cycle
         if (iabs(itype(i+1)).eq.20) iblock=2
         if (iabs(itype(i+1)).ne.20) iblock=1
         dethetai=0.0d0
@@ -3436,7 +3447,7 @@ c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
           coskt(k)=dcos(k*theti2)
           sinkt(k)=dsin(k*theti2)
         enddo
-        if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
+        if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
 #ifdef OSF
           phii=phi(i)
           if (phii.ne.phii) phii=150.0
@@ -3450,13 +3461,14 @@ c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
           enddo
         else
           phii=0.0d0
-          ityp1=nthetyp+1
+c          ityp1=nthetyp+1
           do k=1,nsingle
+            ityp1=ithetyp((itype(i-2)))
             cosph1(k)=0.0d0
             sinph1(k)=0.0d0
           enddo 
         endif
-        if (i.lt.nres .and. itype(i).ne.ntyp1) then
+        if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
 #ifdef OSF
           phii1=phi(i+1)
           if (phii1.ne.phii1) phii1=150.0
@@ -3471,7 +3483,8 @@ c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
           enddo
         else
           phii1=0.0d0
-          ityp3=nthetyp+1
+c          ityp3=nthetyp+1
+          ityp3=ithetyp((itype(i)))
           do k=1,nsingle
             cosph2(k)=0.0d0
             sinph2(k)=0.0d0
@@ -3588,7 +3601,8 @@ c        call flush(iout)
         etheta=etheta+ethetai
         if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii
         if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1
-        gloc(nphi+i-2,icg)=wang*dethetai
+c        gloc(nphi+i-2,icg)=wang*dethetai
+        gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai
       enddo
       return
       end
@@ -3939,7 +3953,7 @@ C     &   dc_norm(3,i+nres)
           y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
         enddo
         do j = 1,3
-          z_prime(j) = -uz(j,i-1)*dsign(1.0,dfloat(itype(i)))
+          z_prime(j) = -uz(j,i-1)*dsign(1.0d0,dfloat(itype(i)))
         enddo     
 c       write (2,*) "i",i
 c       write (2,*) "x_prime",(x_prime(j),j=1,3)
@@ -3979,7 +3993,7 @@ c        write (2,*) "xx",xx," yy",yy," zz",zz
 Cc diagnostics - remove later
         xx1 = dcos(alph(2))
         yy1 = dsin(alph(2))*dcos(omeg(2))
-        zz1 = -dsign(1.0,itype(i))*dsin(alph(2))*dsin(omeg(2))
+        zz1 = -dsign(1.0d0,itype(i))*dsin(alph(2))*dsin(omeg(2))
         write(2,'(3f8.1,3f9.3,1x,3f9.3)') 
      &    alph(2)*rad2deg,omeg(2)*rad2deg,theta(3)*rad2deg,xx,yy,zz,
      &    xx1,yy1,zz1
@@ -4022,6 +4036,8 @@ c     &   dscp1,dscp2,sumene
 c        sumene = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
         escloc = escloc + sumene
 c        write (2,*) "escloc",escloc
+c        write (2,*) "i",i," escloc",sumene,escloc,it,itype(i),
+c     &  zz,xx,yy
         if (.not. calc_grad) goto 1
 #ifdef DEBUG
 C
@@ -4151,9 +4167,9 @@ c     &   (dC_norm(j,i-1),j=1,3)," vbld_inv",vbld_inv(i+1),vbld_inv(i)
          dZZ_Ci(k)=0.0d0
          do j=1,3
            dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1)
-     & *dsign(1.0,dfloat(itype(i)))*dC_norm(j,i+nres)
+     & *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres)
            dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1)
-     &  *dsign(1.0,dfloat(itype(i)))*dC_norm(j,i+nres)
+     &  *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres)
          enddo
           
          dXX_XYZ(k)=vbld_inv(i+nres)*(x_prime(k)-xx*dC_norm(k,i+nres))
@@ -4724,9 +4740,9 @@ c------------------------------------------------------------------------------
       integer dimen1,dimen2,atom,indx
       double precision buffer(dimen1,dimen2)
       double precision zapas 
-      common /contacts_hb/ zapas(3,20,maxres,7),
-     &   facont_hb(20,maxres),ees0p(20,maxres),ees0m(20,maxres),
-     &         num_cont_hb(maxres),jcont_hb(20,maxres)
+      common /contacts_hb/ zapas(3,ntyp,maxres,7),
+     &   facont_hb(ntyp,maxres),ees0p(ntyp,maxres),ees0m(ntyp,maxres),
+     &         num_cont_hb(maxres),jcont_hb(ntyp,maxres)
       num_kont=num_cont_hb(atom)
       do i=1,num_kont
         do k=1,7
@@ -4749,9 +4765,10 @@ c------------------------------------------------------------------------------
       integer dimen1,dimen2,atom,indx
       double precision buffer(dimen1,dimen2)
       double precision zapas 
-      common /contacts_hb/ zapas(3,20,maxres,7),
-     &         facont_hb(20,maxres),ees0p(20,maxres),ees0m(20,maxres),
-     &         num_cont_hb(maxres),jcont_hb(20,maxres)
+      common /contacts_hb/ zapas(3,ntyp,maxres,7),
+     &         facont_hb(ntyp,maxres),ees0p(ntyp,maxres),
+     &         ees0m(ntyp,maxres),
+     &         num_cont_hb(maxres),jcont_hb(ntyp,maxres)
       num_kont=buffer(1,indx+26)
       num_kont_old=num_cont_hb(atom)
       num_cont_hb(atom)=num_kont+num_kont_old
@@ -6490,7 +6507,7 @@ c----------------------------------------------------------------------------
       include 'COMMON.GEO'
       logical swap
       double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2),
-     & auxvec1(2),auxvec2(1),auxmat1(2,2)
+     & auxvec1(2),auxvec2(2),auxmat1(2,2)
       logical lprn
       common /kutas/ lprn
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC