projects / unres.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
Dostosowanie w src whama i clustra
[unres.git] / source / wham / src / energy_p_new.F
diff --git a/source/wham/src/energy_p_new.F b/source/wham/src/energy_p_new.F
index 8303bc2..75d67e1 100644 (file)
--- a/source/wham/src/energy_p_new.F
+++ b/source/wham/src/energy_p_new.F
@@ -228,7 +228,7 @@ C
      &   +wturn3*fact(2)*gel_loc_turn3(i)
      &   +wturn6*fact(5)*gel_loc_turn6(i)
      &   +wel_loc*fact(2)*gel_loc_loc(i)
-     &   +wsccor*fact(1)*gsccor_loc(i)
+c     &   +wsccor*fact(1)*gsccor_loc(i)
       enddo
       endif
       return
@@ -369,6 +369,7 @@ cd    print *,'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
       evdw_t=0.0d0
       do i=iatsc_s,iatsc_e
         itypi=iabs(itype(i))
+        if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
         xi=c(1,nres+i)
         yi=c(2,nres+i)
@@ -383,6 +384,7 @@ cd        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
 cd   &                  'iend=',iend(i,iint)
           do j=istart(i,iint),iend(i,iint)
             itypj=iabs(itype(j))
+            if (itypj.eq.ntyp1) cycle
             xj=c(1,nres+j)-xi
             yj=c(2,nres+j)-yi
             zj=c(3,nres+j)-zi
@@ -540,6 +542,7 @@ c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
       evdw_t=0.0d0
       do i=iatsc_s,iatsc_e
         itypi=iabs(itype(i))
+        if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
         xi=c(1,nres+i)
         yi=c(2,nres+i)
@@ -550,6 +553,7 @@ C
         do iint=1,nint_gr(i)
           do j=istart(i,iint),iend(i,iint)
             itypj=iabs(itype(j))
+            if (itypj.eq.ntyp1) cycle
             xj=c(1,nres+j)-xi
             yj=c(2,nres+j)-yi
             zj=c(3,nres+j)-zi
@@ -3436,11 +3440,13 @@ C
       etheta=0.0D0
 c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
       do i=ithet_start,ithet_end
+        if (iabs(itype(i+1)).eq.20) iblock=2
+        if (iabs(itype(i+1)).ne.20) iblock=1
         dethetai=0.0d0
         dephii=0.0d0
         dephii1=0.0d0
         theti2=0.5d0*theta(i)
-        ityp2=ithetyp(iabs(itype(i-1)))
+        ityp2=ithetyp((itype(i-1)))
         do k=1,nntheterm
           coskt(k)=dcos(k*theti2)
           sinkt(k)=dsin(k*theti2)
@@ -3473,7 +3479,7 @@ c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
 #else
           phii1=phi(i+1)
 #endif
-          ityp3=ithetyp(iabs(itype(i)))
+          ityp3=ithetyp((itype(i)))
           do k=1,nsingle
             cosph2(k)=dcos(k*phii1)
             sinph2(k)=dsin(k*phii1)
@@ -3489,7 +3495,7 @@ c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
 c        write (iout,*) "i",i," ityp1",itype(i-2),ityp1,
 c     &   " ityp2",itype(i-1),ityp2," ityp3",itype(i),ityp3
 c        call flush(iout)
-        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)
@@ -3511,11 +3517,12 @@ c        call flush(iout)
         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
@@ -3534,24 +3541,24 @@ c        call flush(iout)
         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)
@@ -3559,28 +3566,29 @@ c        call flush(iout)
         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)
@@ -3944,7 +3952,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)
@@ -3966,7 +3974,7 @@ c
         do j = 1,3
           xx = xx + x_prime(j)*dc_norm(j,i+nres)
           yy = yy + y_prime(j)*dc_norm(j,i+nres)
-          zz = zz + dsign(1.0,itype(i))*z_prime(j)*dc_norm(j,i+nres)
+          zz = zz + z_prime(j)*dc_norm(j,i+nres)
         enddo
 
         xxtab(i)=xx
@@ -3984,7 +3992,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
@@ -4155,8 +4163,11 @@ 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))
@@ -4621,7 +4632,7 @@ c     &gloc_sc(intertyp,i-3,icg)
      &  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)
-        gsccor_loc(i-3)=gsccor_loc(i-3)+gloci
+c        gsccor_loc(i-3)=gsccor_loc(i-3)+gloci
        enddo !intertyp
       enddo
 c        do i=1,nres
@@ -4734,9 +4745,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
@@ -4759,9 +4770,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
UNRESPACK 3.x