By Adam
[unres.git] / source / unres / src_MD / intcartderiv.F
index d029e18..889bedd 100644 (file)
@@ -12,6 +12,7 @@
       include 'COMMON.DERIV'
       include 'COMMON.IOUNITS'
       include 'COMMON.LOCAL'
+      include 'COMMON.SCCOR'
       double precision dcostheta(3,2,maxres),
      & dcosphi(3,3,maxres),dsinphi(3,3,maxres),
      & dcosalpha(3,3,maxres),dcosomega(3,3,maxres),
@@ -51,20 +52,20 @@ c We need dtheta(:,:,i-1) to compute dphi(:,:,i)
 #else
       do i=3,nres
 #endif
-      if (itype(i-1).ne.10) then
+      if ((itype(i-1).ne.10).and.(itype(i-1).ne.21)) then
         cost1=dcos(omicron(1,i))
        sint1=sqrt(1-cost1*cost1)
         cost2=dcos(omicron(2,i))
         sint2=sqrt(1-cost2*cost2)
         do j=1,3
 CC Calculate derivative over first omicron (Cai-2,Cai-1,SCi-1) 
-          dcosomicron(j,1,1,i)=-(-dc_norm(j,i-1+nres)+
+          dcosomicron(j,1,1,i)=-(dc_norm(j,i-1+nres)+
      &    cost1*dc_norm(j,i-2))/
      &   vbld(i-1)
           domicron(j,1,1,i)=-1/sint1*dcosomicron(j,1,1,i)     
           dcosomicron(j,1,2,i)=-(dc_norm(j,i-2)
-     &    +cost1*(-dc_norm(j,i-1+nres)))/
-     &   vbld(i+nres)
+     &    +cost1*(dc_norm(j,i-1+nres)))/
+     &   vbld(i-1+nres)
           domicron(j,1,2,i)=-1/sint1*dcosomicron(j,1,2,i)  
 CC Calculate derivative over second omicron Sci-1,Cai-1 Cai
 CC Looks messy but better than if in loop
@@ -73,8 +74,9 @@ CC Looks messy but better than if in loop
      &    vbld(i)
           domicron(j,2,1,i)=-1/sint2*dcosomicron(j,2,1,i)
           dcosomicron(j,2,2,i)=-(dc_norm(j,i-1)
-     &     +cost1*(-dc_norm(j,i-1+nres)))/
-     &    vbld(i+nres)
+     &     +cost2*(-dc_norm(j,i-1+nres)))/
+     &    vbld(i-1+nres)
+c          write(iout,*) "vbld", i,itype(i),vbld(i-1+nres)
           domicron(j,2,2,i)=-1/sint2*dcosomicron(j,2,2,i)   
         enddo
        endif
@@ -147,12 +149,18 @@ c   Obtaining the gamma derivatives from cosine derivative
         endif                                                                                           
       enddo
 
+      do i=1,nres-1
+       do j=1,3
+        dc_norm2(j,i+nres)=-dc_norm(j,i+nres)
+       enddo
+      enddo
 Calculate derivative of Tauangle
 #ifdef PARINTDER
-      do i=iphi1_start-1,iphi1_end
+      do i=itau_start,itau_end
 #else
       do i=3,nres
 #endif
+       if ((itype(i-2).eq.21).or.(itype(i-2).eq.10)) cycle
 cc dtauangle(j,intertyp,dervityp,residue number)
 cc INTERTYP=1 SC...Ca...Ca..Ca
 c the conventional case
@@ -162,22 +170,21 @@ c the conventional case
         cost=dcos(theta(i))
         cost1=dcos(omicron(2,i-1))
         cosg=dcos(tauangle(1,i))
-        do j=1,3
-        dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres)
-        enddo
+C        do j=1,3
+C        dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres)
+cc       write(iout,*) dc_norm2(j,i-2+nres),"dcnorm"
+C        enddo
         scalp=scalar(dc_norm2(1,i-2+nres),dc_norm(1,i-1))
         fac0=1.0d0/(sint1*sint)
         fac1=cost*fac0
         fac2=cost1*fac0
         fac3=cosg*cost1/(sint1*sint1)
         fac4=cosg*cost/(sint*sint)
+cc         write(iout,*) "faki",fac0,fac1,fac2,fac3,fac4
 c    Obtaining the gamma derivatives from sine derivative                                
        if (tauangle(1,i).gt.-pi4.and.tauangle(1,i).le.pi4.or.
      &     tauangle(1,i).gt.pi34.and.tauangle(1,i).le.pi.or.
      &     tauangle(1,i).gt.-pi.and.tauangle(1,i).le.-pi34) then
-         do j=1,3
-         dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres)
-         enddo
          call vecpr(dc_norm(1,i-1),dc_norm(1,i-2),vp1)
          call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-1),vp2)
          call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-2),vp3)
@@ -185,55 +192,61 @@ c    Obtaining the gamma derivatives from sine derivative
             ctgt=cost/sint
             ctgt1=cost1/sint1
             cosg_inv=1.0d0/cosg
-            dsintau(j,1,1,i)=-sing*ctgt1*domicron(j,2,1,i-1)
-     &-(fac0*vp1(j)+sing*(-dc_norm(j,i-2+nres)))
+            dsintau(j,1,1,i)=-sing*ctgt1*domicron(j,2,2,i-1)
+     &-(fac0*vp1(j)+sing*(dc_norm2(j,i-2+nres)))
      & *vbld_inv(i-2+nres)
             dtauangle(j,1,1,i)=cosg_inv*dsintau(j,1,1,i)
             dsintau(j,1,2,i)=
-     &        -sing*(ctgt1*domicron(j,2,2,i-1)+ctgt*dtheta(j,1,i))
+     &        -sing*(ctgt1*domicron(j,2,1,i-1)+ctgt*dtheta(j,1,i))
      &        -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
+c            write(iout,*) "dsintau", dsintau(j,1,2,i)
             dtauangle(j,1,2,i)=cosg_inv*dsintau(j,1,2,i)
 c Bug fixed 3/24/05 (AL)
             dsintau(j,1,3,i)=-sing*ctgt*dtheta(j,2,i)
      &        +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i)
 c     &        +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
-            dtauangle(j,1,3,i)=cosg_inv*dsinphi(j,1,3,i)
+            dtauangle(j,1,3,i)=cosg_inv*dsintau(j,1,3,i)
          enddo                                          
 c   Obtaining the gamma derivatives from cosine derivative
         else
            do j=1,3
-           dcostau(j,1,1,i)=fac1*dcosomicron(j,2,1,i-1)+fac3*
-     &     dcosomicron(j,2,1,i-1)-fac0*(dc_norm(j,i-1)-scalp*
-     &     (-dc_norm(j,i-2+nres)))/vbld(i-2+nres)
+           dcostau(j,1,1,i)=fac1*dcosomicron(j,2,2,i-1)+fac3*
+     &     dcosomicron(j,2,2,i-1)-fac0*(dc_norm(j,i-1)-scalp*
+     &     (dc_norm2(j,i-2+nres)))/vbld(i-2+nres)
            dtauangle(j,1,1,i)=-1/sing*dcostau(j,1,1,i)
-           dcostau(j,1,2,i)=fac1*dcosomicron(j,2,2,i-1)+fac2*
-     &     dcostheta(j,1,i)+fac3*dcosomicron(j,2,2,i-1)+fac4*
+           dcostau(j,1,2,i)=fac1*dcosomicron(j,2,1,i-1)+fac2*
+     &     dcostheta(j,1,i)+fac3*dcosomicron(j,2,1,i-1)+fac4*
      &     dcostheta(j,1,i)
            dtauangle(j,1,2,i)=-1/sing*dcostau(j,1,2,i)
            dcostau(j,1,3,i)=fac2*dcostheta(j,2,i)+fac4*
      &     dcostheta(j,2,i)-fac0*(-dc_norm(j,i-2+nres)-scalp*
      &     dc_norm(j,i-1))/vbld(i)
            dtauangle(j,1,3,i)=-1/sing*dcostau(j,1,3,i)
+c         write (iout,*) "else",i
          enddo
-        endif                                                                                            
+        endif
+c        do k=1,3                 
+c        write(iout,*) "tu",i,k,(dtauangle(j,1,k,i),j=1,3)        
+c        enddo                
       enddo
 CC Second case Ca...Ca...Ca...SC
 #ifdef PARINTDER
-      do i=iphi1_start,iphi1_end
+      do i=itau_start,itau_end
 #else
       do i=4,nres
 #endif
+       if ((itype(i-1).eq.21).or.(itype(i-1).eq.10)) cycle
 c the conventional case
         sint=dsin(omicron(1,i))
         sint1=dsin(theta(i-1))
-        sing=dsin(tauangle(2,i)
-        cost=dcos(omicron(1,i)
+        sing=dsin(tauangle(2,i))
+        cost=dcos(omicron(1,i))
         cost1=dcos(theta(i-1))
         cosg=dcos(tauangle(2,i))
-        do j=1,3
-        dc_norm2(j,i-1+nres)=-dc_norm(j,i-1+nres)
-        enddo
-        scalp=scalar(dc_norm(1,i-3),dc_norm2(1,i-1+nres))
+c        do j=1,3
+c        dc_norm2(j,i-1+nres)=-dc_norm(j,i-1+nres)
+c        enddo
+        scalp=scalar(dc_norm(1,i-3),dc_norm(1,i-1+nres))
         fac0=1.0d0/(sint1*sint)
         fac1=cost*fac0
         fac2=cost1*fac0
@@ -244,40 +257,46 @@ c    Obtaining the gamma derivatives from sine derivative
      &     tauangle(2,i).gt.pi34.and.tauangle(2,i).le.pi.or.
      &     tauangle(2,i).gt.-pi.and.tauangle(2,i).le.-pi34) then
          call vecpr(dc_norm2(1,i-1+nres),dc_norm(1,i-2),vp1)
-         call vecpr(dc_norm(1,i-3),dc_norm2(1,i-1+nres),vp2)
+         call vecpr(dc_norm(1,i-3),dc_norm(1,i-1+nres),vp2)
          call vecpr(dc_norm(1,i-3),dc_norm(1,i-2),vp3)
         do j=1,3
             ctgt=cost/sint
             ctgt1=cost1/sint1
             cosg_inv=1.0d0/cosg
             dsintau(j,2,1,i)=-sing*ctgt1*dtheta(j,1,i-1)
-     &        -(fac0*vp1(j)+sing*dc_norm(j,i-3))*vbld_inv(i-2)
+     &        +(fac0*vp1(j)-sing*dc_norm(j,i-3))*vbld_inv(i-2)
+c       write(iout,*) i,j,dsintau(j,2,1,i),sing*ctgt1*dtheta(j,1,i-1),
+c     &fac0*vp1(j),sing*dc_norm(j,i-3),vbld_inv(i-2),"dsintau(2,1)"
             dtauangle(j,2,1,i)=cosg_inv*dsintau(j,2,1,i)
             dsintau(j,2,2,i)=
      &        -sing*(ctgt1*dtheta(j,2,i-1)+ctgt*domicron(j,1,1,i))
      &        -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
-            dphi(j,2,i)=cosg_inv*dsintau(j,2,2,i)
+c            write(iout,*) "sprawdzenie",i,j,sing*ctgt1*dtheta(j,2,i-1),
+c     & sing*ctgt*domicron(j,1,2,i),
+c     & (fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
+            dtauangle(j,2,2,i)=cosg_inv*dsintau(j,2,2,i)
 c Bug fixed 3/24/05 (AL)
             dsintau(j,2,3,i)=-sing*ctgt*domicron(j,1,2,i)
-     &        +(fac0*vp3(j)-sing*dc_norm2(j,i-1+nres))*vbld_inv(i)
+     &       +(fac0*vp3(j)-sing*dc_norm(j,i-1+nres))*vbld_inv(i-1+nres)
 c     &        +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
-            dtauangle(j,2,3,i)=cosg_inv*dsinphi(j,2,3,i)
+            dtauangle(j,2,3,i)=cosg_inv*dsintau(j,2,3,i)
          enddo                                          
 c   Obtaining the gamma derivatives from cosine derivative
         else
            do j=1,3
            dcostau(j,2,1,i)=fac1*dcostheta(j,1,i-1)+fac3*
-     &     dcostheta(j,1,i-1)-fac0*(dc_norm2(j,i-1+nres)-scalp*
+     &     dcostheta(j,1,i-1)-fac0*(dc_norm(j,i-1+nres)-scalp*
      &     dc_norm(j,i-3))/vbld(i-2)
            dtauangle(j,2,1,i)=-1/sing*dcostau(j,2,1,i)
            dcostau(j,2,2,i)=fac1*dcostheta(j,2,i-1)+fac2*
      &     dcosomicron(j,1,1,i)+fac3*dcostheta(j,2,i-1)+fac4*
      &     dcosomicron(j,1,1,i)
-           dtauanlge(j,2,2,i)=-1/sing*dcostau(j,2,2,i)
+           dtauangle(j,2,2,i)=-1/sing*dcostau(j,2,2,i)
            dcostau(j,2,3,i)=fac2*dcosomicron(j,1,2,i)+fac4*
-     &     dcostheta(j,1,2,i)-fac0*(dc_norm(j,i-3)-scalp*
-     &     dc_norm2(j,i-1+nres))/vbld(i-1+nres)
-           dtauanlge(j,2,3,i)=-1/sing*dcosphi(j,3,i)
+     &     dcosomicron(j,1,2,i)-fac0*(dc_norm(j,i-3)-scalp*
+     &     dc_norm(j,i-1+nres))/vbld(i-1+nres)
+           dtauangle(j,2,3,i)=-1/sing*dcostau(j,2,3,i)
+c        write(iout,*) i,j,"else", dtauangle(j,2,3,i) 
          enddo
         endif                                                                                            
       enddo
@@ -286,22 +305,24 @@ c   Obtaining the gamma derivatives from cosine derivative
 CCC third case SC...Ca...Ca...SC
 #ifdef PARINTDER
 
-      do i=iphi1_start-1,iphi1_end
+      do i=itau_start,itau_end
 #else
       do i=3,nres
 #endif
 c the conventional case
+      if ((itype(i-1).eq.21).or.(itype(i-1).eq.10).or.
+     &(itype(i-2).eq.21).or.(itype(i-2).eq.10)) cycle
         sint=dsin(omicron(1,i))
         sint1=dsin(omicron(2,i-1))
         sing=dsin(tauangle(3,i))
         cost=dcos(omicron(1,i))
         cost1=dcos(omicron(2,i-1))
         cosg=dcos(tauangle(3,i))
-        do j=1,3
-        dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres)
-        dc_norm2(j,i-1+nres)=-dc_norm(j,i-2+nres)
-        enddo
-        scalp=scalar(dc_norm2(1,i-2+nres),dc_norm2(1,i-1+nres))
+C        do j=1,3
+C        dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres)
+c        dc_norm2(j,i-1+nres)=-dc_norm(j,i-1+nres)
+C        enddo
+        scalp=scalar(dc_norm2(1,i-2+nres),dc_norm(1,i-1+nres))
         fac0=1.0d0/(sint1*sint)
         fac1=cost*fac0
         fac2=cost1*fac0
@@ -311,43 +332,44 @@ c    Obtaining the gamma derivatives from sine derivative
        if (tauangle(3,i).gt.-pi4.and.tauangle(3,i).le.pi4.or.
      &     tauangle(3,i).gt.pi34.and.tauangle(3,i).le.pi.or.
      &     tauangle(3,i).gt.-pi.and.tauangle(3,i).le.-pi34) then
-         call vecpr(dc_norm2(1,i-1+nres),dc_norm(1,i-2),vp1)
-         call vecpr(dc_norm2(1,i-2+nres),dc_norm2(1,i-1+nres),vp2)
+         call vecpr(dc_norm(1,i-1+nres),dc_norm(1,i-2),vp1)
+         call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-1+nres),vp2)
          call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-2),vp3)
         do j=1,3
             ctgt=cost/sint
             ctgt1=cost1/sint1
             cosg_inv=1.0d0/cosg
-            dsintau(j,3,1,i)=-sing*ctgt1*domicron(j,2,1,i-1)
-     &        -(fac0*vp1(j)+sing*dc_norm2(j,i-2+nres))
+            dsintau(j,3,1,i)=-sing*ctgt1*domicron(j,2,2,i-1)
+     &        -(fac0*vp1(j)-sing*dc_norm(j,i-2+nres))
      &        *vbld_inv(i-2+nres)
             dtauangle(j,3,1,i)=cosg_inv*dsintau(j,3,1,i)
             dsintau(j,3,2,i)=
-     &        -sing*(ctgt1*domicron(j,2,2,i-1)+ctgt*domicron(j,1,1,i))
-     &        -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1+nres)
+     &        -sing*(ctgt1*domicron(j,2,1,i-1)+ctgt*domicron(j,1,1,i))
+     &        -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
             dtauangle(j,3,2,i)=cosg_inv*dsintau(j,3,2,i)
 c Bug fixed 3/24/05 (AL)
             dsintau(j,3,3,i)=-sing*ctgt*domicron(j,1,2,i)
-     &        +(fac0*vp3(j)-sing*dc_norm2(j,i-1+nres))
+     &        +(fac0*vp3(j)-sing*dc_norm(j,i-1+nres))
      &        *vbld_inv(i-1+nres)
 c     &        +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
-            dphi(j,3,3,i)=cosg_inv*dsintau(j,3,3,i)
+            dtauangle(j,3,3,i)=cosg_inv*dsintau(j,3,3,i)
          enddo                                          
 c   Obtaining the gamma derivatives from cosine derivative
         else
            do j=1,3
-           dcostau(j,3,1,i)=fac1*dcosomicron(j,2,1,i-1)+fac3*
-     &     dcostheta(j,1,i-1)-fac0*(dc_norm2(j,i-1+nres)-scalp*
+           dcostau(j,3,1,i)=fac1*dcosomicron(j,2,2,i-1)+fac3*
+     &     dcosomicron(j,2,2,i-1)-fac0*(dc_norm(j,i-1+nres)-scalp*
      &     dc_norm2(j,i-2+nres))/vbld(i-2+nres)
            dtauangle(j,3,1,i)=-1/sing*dcostau(j,3,1,i)
-           dcostau(j,3,2,i)=fac1*dcosomicron(j,2,2,i-1)+fac2*
-     &     dcosomicron(j,1,1,i)+fac3*dcosomicron(j,2,2,i-1)+fac4*
+           dcostau(j,3,2,i)=fac1*dcosomicron(j,2,1,i-1)+fac2*
+     &     dcosomicron(j,1,1,i)+fac3*dcosomicron(j,2,1,i-1)+fac4*
      &     dcosomicron(j,1,1,i)
            dtauangle(j,3,2,i)=-1/sing*dcostau(j,3,2,i)
            dcostau(j,3,3,i)=fac2*dcosomicron(j,1,2,i)+fac4*
-     &     dcostau(j,3,2,i)-fac0*(dc_norm2(j,i-2+nres)-scalp*
-     &     dc_norm2(j,i-1+nres))/vbld(i-1+nres)
+     &     dcosomicron(j,1,2,i)-fac0*(dc_norm2(j,i-2+nres)-scalp*
+     &     dc_norm(j,i-1+nres))/vbld(i-1+nres)
            dtauangle(j,3,3,i)=-1/sing*dcostau(j,3,3,i)
+c          write(iout,*) "else",i 
          enddo
         endif                                                                                            
       enddo
@@ -441,12 +463,12 @@ c   obtaining the derivatives of omega from cosines
 #if defined(MPI) && defined(PARINTDER)
       if (nfgtasks.gt.1) then
 #ifdef DEBUG
-cd      write (iout,*) "Gather dtheta"
+       write (iout,*) "Gather dtheta"
 cd      call flush(iout)
-      write (iout,*) "dtheta before gather"
-      do i=1,nres
-        write (iout,'(i3,3(3f8.5,3x))') i,((dtheta(j,k,i),k=1,3),j=1,2)
-      enddo
+c      write (iout,*) "dtheta before gather"
+c      do i=1,nres
+c        write (iout,'(i3,3(3f8.5,3x))') i,((dtheta(j,k,i),k=1,3),j=1,2)
+c      enddo
 #endif
       call MPI_Gatherv(dtheta(1,1,ithet_start),ithet_count(fg_rank),
      &  MPI_THET,dtheta(1,1,1),ithet_count(0),ithet_displ(0),MPI_THET,