pierwsza proba merga devela do adasko
[unres.git] / source / unres / src_MD / energy_p_new_barrier.F
index a2eb0f9..74bd3ae 100644 (file)
@@ -1449,7 +1449,7 @@ C
         do iint=1,nint_gr(i)
           do j=istart(i,iint),iend(i,iint)
             ind=ind+1
-            itypj=itype(j)
+            itypj=iabs(itype(j))
 c            dscj_inv=dsc_inv(itypj)
             dscj_inv=vbld_inv(j+nres)
             chi1=chi(itypi,itypj)
@@ -4808,6 +4808,9 @@ C
           enddo 
         endif
         if (i.lt.nres) then
+
+        if (iabs(itype(i+1)).eq.20) iblock=2
+        if (iabs(itype(i+1)).ne.20) iblock=1
 #ifdef OSF
           phii1=phi(i+1)
           if (phii1.ne.phii1) phii1=150.0
@@ -4828,7 +4831,7 @@ C
             sinph2(k)=0.0d0
           enddo
         endif  
-        ethetai=aa0thet(ityp1,ityp2,ityp3)
+         ethetai=aa0thet(ityp1,ityp2,ityp3,iblock)
         do k=1,ndouble
           do l=1,k-1
             ccl=cosph1(l)*cosph2(k-l)
         enddo
         endif
         do k=1,ntheterm
-          ethetai=ethetai+aathet(k,ityp1,ityp2,ityp3)*sinkt(k)
-          dethetai=dethetai+0.5d0*k*aathet(k,ityp1,ityp2,ityp3)
+          ethetai=ethetai+aathet(k,ityp1,ityp2,ityp3,iblock)*sinkt(k)
+          dethetai=dethetai+0.5d0*k*aathet(k,ityp1,ityp2,ityp3,iblock)
      &      *coskt(k)
           if (lprn)
-     &    write (iout,*) "k",k," aathet",aathet(k,ityp1,ityp2,ityp3),
+     &    write (iout,*) "k",k,
+     &    "aathet",aathet(k,ityp1,ityp2,ityp3,iblock),
      &     " ethetai",ethetai
         enddo
         if (lprn) then
         endif
         do m=1,ntheterm2
           do k=1,nsingle
-            aux=bbthet(k,m,ityp1,ityp2,ityp3)*cosph1(k)
-     &         +ccthet(k,m,ityp1,ityp2,ityp3)*sinph1(k)
-     &         +ddthet(k,m,ityp1,ityp2,ityp3)*cosph2(k)
-     &         +eethet(k,m,ityp1,ityp2,ityp3)*sinph2(k)
+            aux=bbthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph1(k)
+     &         +ccthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph1(k)
+     &         +ddthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph2(k)
+     &         +eethet(k,m,ityp1,ityp2,ityp3,iblock)*sinph2(k)
             ethetai=ethetai+sinkt(m)*aux
             dethetai=dethetai+0.5d0*m*aux*coskt(m)
             dephii=dephii+k*sinkt(m)*(
-     &          ccthet(k,m,ityp1,ityp2,ityp3)*cosph1(k)-
-     &          bbthet(k,m,ityp1,ityp2,ityp3)*sinph1(k))
+     &          ccthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph1(k)-
+     &          bbthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph1(k))
             dephii1=dephii1+k*sinkt(m)*(
-     &          eethet(k,m,ityp1,ityp2,ityp3)*cosph2(k)-
-     &          ddthet(k,m,ityp1,ityp2,ityp3)*sinph2(k))
+     &          eethet(k,m,ityp1,ityp2,ityp3,iblock)*cosph2(k)-
+     &          ddthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph2(k))
             if (lprn)
      &      write (iout,*) "m",m," k",k," bbthet",
-     &         bbthet(k,m,ityp1,ityp2,ityp3)," ccthet",
-     &         ccthet(k,m,ityp1,ityp2,ityp3)," ddthet",
-     &         ddthet(k,m,ityp1,ityp2,ityp3)," eethet",
-     &         eethet(k,m,ityp1,ityp2,ityp3)," ethetai",ethetai
+     &         bbthet(k,m,ityp1,ityp2,ityp3,iblock)," ccthet",
+     &         ccthet(k,m,ityp1,ityp2,ityp3,iblock)," ddthet",
+     &         ddthet(k,m,ityp1,ityp2,ityp3,iblock)," eethet",
+     &         eethet(k,m,ityp1,ityp2,ityp3,iblock)," ethetai",ethetai
           enddo
         enddo
         if (lprn)
         do m=1,ntheterm3
           do k=2,ndouble
             do l=1,k-1
-              aux=ffthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+
-     &            ffthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)+
-     &            ggthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+
-     &            ggthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)
+       aux=ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)+
+     & ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l)+
+     & ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)+
+     & ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l)
+
               ethetai=ethetai+sinkt(m)*aux
               dethetai=dethetai+0.5d0*m*coskt(m)*aux
               dephii=dephii+l*sinkt(m)*(
-     &           -ffthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)-
-     &            ffthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+
-     &            ggthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+
-     &            ggthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l))
+     & -ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)-
+     &  ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l)+
+     &  ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)+
+     &  ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l))
+
               dephii1=dephii1+(k-l)*sinkt(m)*(
-     &           -ffthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+
-     &            ffthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+
-     &            ggthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)-
-     &            ggthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l))
+     &-ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)+
+     & ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l)+
+     & ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)-
+     & ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l))
+
               if (lprn) then
               write (iout,*) "m",m," k",k," l",l," ffthet",
-     &            ffthet(l,k,m,ityp1,ityp2,ityp3),
-     &            ffthet(k,l,m,ityp1,ityp2,ityp3)," ggthet",
-     &            ggthet(l,k,m,ityp1,ityp2,ityp3),
-     &            ggthet(k,l,m,ityp1,ityp2,ityp3)," ethetai",ethetai
+     &            ffthet(l,k,m,ityp1,ityp2,ityp3,iblock),
+     &            ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)," ggthet",
+     &            ggthet(l,k,m,ityp1,ityp2,ityp3,iblock),
+     &            ggthet(k,l,m,ityp1,ityp2,ityp3,iblock),
+     &            " ethetai",ethetai
+
               write (iout,*) cosph1ph2(l,k)*sinkt(m),
      &            cosph1ph2(k,l)*sinkt(m),
      &            sinph1ph2(l,k)*sinkt(m),sinph1ph2(k,l)*sinkt(m)
@@ -5271,7 +5280,7 @@ C
         cosfac=dsqrt(cosfac2)
         sinfac2=0.5d0/(1.0d0-costtab(i+1))
         sinfac=dsqrt(sinfac2)
-        it=itype(i)
+        it=iabs(itype(i))
         if (it.eq.10) goto 1
 c
 C  Compute the axes of tghe local cartesian coordinates system; store in
@@ -5289,7 +5298,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)
+          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)
@@ -5321,7 +5330,7 @@ C
 C Compute the energy of the ith side cbain
 C
 c        write (2,*) "xx",xx," yy",yy," zz",zz
-        it=itype(i)
+        it=iabs(itype(i))
         do j = 1,65
           x(j) = sc_parmin(j,it) 
         enddo
@@ -5329,7 +5338,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 = -dsin(alph(2))*dsin(omeg(2))
+        zz1 = -dsign(1.0, dfloat(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
@@ -5499,8 +5508,10 @@ c     &   (dC_norm(j,i-1),j=1,3)," vbld_inv",vbld_inv(i+1),vbld_inv(i)
          dZZ_Ci1(k)=0.0d0
          dZZ_Ci(k)=0.0d0
          do j=1,3
-           dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1)*dC_norm(j,i+nres)
-           dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1)*dC_norm(j,i+nres)
+           dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1)
+     &     *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.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))
@@ -5786,6 +5797,8 @@ c     lprn=.true.
       etors=0.0D0
       do i=iphi_start,iphi_end
       etors_ii=0.0D0
+c        if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1
+c     &      .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
         itori=itortyp(itype(i-2))
         itori1=itortyp(itype(i-1))
         phii=phi(i)
@@ -5879,6 +5892,8 @@ C Set lprn=.true. for debugging
 c     lprn=.true.
       etors_d=0.0D0
       do i=iphid_start,iphid_end
+c        if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1
+c     &      .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
         itori=itortyp(itype(i-2))
         itori1=itortyp(itype(i-1))
         itori2=itortyp(itype(i))
@@ -5951,10 +5966,11 @@ c        amino-acid residues.
 C Set lprn=.true. for debugging
       lprn=.false.
 c      lprn=.true.
-c      write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor
+c     write (iout,*) "EBACK_SC_COR",itau_start,itau_end
       esccor=0.0D0
       do i=itau_start,itau_end
         esccor_ii=0.0D0
+        if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle
         isccori=isccortyp(itype(i-2))
         isccori1=isccortyp(itype(i-1))
         phii=phi(i)
@@ -5973,14 +5989,16 @@ c   2 = Ca...Ca...Ca...SC
 c   3 = SC...Ca...Ca...SCi
         gloci=0.0D0
         if (((intertyp.eq.3).and.((itype(i-2).eq.10).or.
-     &      (itype(i-1).eq.10).or.(itype(i-2).eq.21).or.
-     &      (itype(i-1).eq.21)))
+     &      (itype(i-1).eq.10).or.(itype(i-2).eq.ntyp1).or.
+     &      (itype(i-1).eq.ntyp1)))
      &    .or. ((intertyp.eq.1).and.((itype(i-2).eq.10)
-     &     .or.(itype(i-2).eq.21)))
+     &     .or.(itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)
+     &     .or.(itype(i).eq.ntyp1)))
      &    .or.((intertyp.eq.2).and.((itype(i-1).eq.10).or.
-     &      (itype(i-1).eq.21)))) cycle  
-        if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1).eq.21)) cycle
-        if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres).eq.21))
+     &      (itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or.
+     &      (itype(i-3).eq.ntyp1)))) cycle
+        if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1).eq.ntyp1)) cycle
+        if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres).eq.ntyp1))
      & cycle
         do j=1,nterm_sccor(isccori,isccori1)
           v1ij=v1sccor(j,intertyp,isccori,isccori1)
@@ -5995,9 +6013,9 @@ c        write (iout,*) "WTF",intertyp,i,itype(i),v1ij*cosphi+v2ij*sinphi
 c     &gloc_sc(intertyp,i-3,icg)
         if (lprn)
      &  write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)')
-     &  restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,
-     &  (v1sccor(j,intertyp,itori,itori1),j=1,6)
-     & ,(v2sccor(j,intertyp,itori,itori1),j=1,6)
+     &  restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,isccori,isccori1,
+     &  (v1sccor(j,intertyp,isccori,isccori1),j=1,6)
+     & ,(v2sccor(j,intertyp,isccori,isccori1),j=1,6)
         gsccor_loc(i-3)=gsccor_loc(i-3)+gloci
        enddo !intertyp
       enddo