update
[unres.git] / source / unres / src_MD / intcartderiv.F
index c220540..d0c1e84 100644 (file)
@@ -35,7 +35,7 @@ c We need dtheta(:,:,i-1) to compute dphi(:,:,i)
       do i=3,nres
 #endif
         cost=dcos(theta(i))
-       sint=sqrt(1-cost*cost)
+       sint=dsqrt(1-cost*cost)
         do j=1,3
           dcostheta(j,1,i)=-(dc_norm(j,i-1)+cost*dc_norm(j,i-2))/
      &   vbld(i-1)
@@ -54,9 +54,9 @@ c We need dtheta(:,:,i-1) to compute dphi(:,:,i)
 #endif
       if ((itype(i-1).ne.10).and.(itype(i-1).ne.21)) then
         cost1=dcos(omicron(1,i))
-       sint1=sqrt(1-cost1*cost1)
+       sint1=dsqrt(1-cost1*cost1)
         cost2=dcos(omicron(2,i))
-        sint2=sqrt(1-cost2*cost2)
+        sint2=dsqrt(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)+
@@ -109,7 +109,7 @@ c the conventional case
 c    Obtaining the gamma derivatives from sine derivative                               
        if (phi(i).gt.-pi4.and.phi(i).le.pi4.or.
      &     phi(i).gt.pi34.and.phi(i).le.pi.or.
-     &     phi(i).gt.-pi.and.phi(i).le.-pi34) then
+     &     phi(i).ge.-pi.and.phi(i).le.-pi34) then
          call vecpr(dc_norm(1,i-1),dc_norm(1,i-2),vp1)
          call vecpr(dc_norm(1,i-3),dc_norm(1,i-1),vp2)
          call vecpr(dc_norm(1,i-3),dc_norm(1,i-2),vp3) 
@@ -149,6 +149,11 @@ c   Obtaining the gamma derivatives from cosine derivative
         endif                                                                                           
       enddo
 
+      do i=1,nres
+       do j=1,3
+        dc_norm2(j,i+nres)=-dc_norm(j,i+nres)
+       enddo
+      enddo
 Calculate derivative of Tauangle
 #ifdef PARINTDER
       do i=itau_start,itau_end
@@ -165,10 +170,10 @@ 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)
+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"
-        enddo
+C        enddo
         scalp=scalar(dc_norm2(1,i-2+nres),dc_norm(1,i-1))
         fac0=1.0d0/(sint1*sint)
         fac1=cost*fac0
@@ -179,7 +184,7 @@ 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
+     &     tauangle(1,i).ge.-pi.and.tauangle(1,i).le.-pi34) then
          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)
@@ -208,15 +213,15 @@ c   Obtaining the gamma derivatives from cosine derivative
            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)
+           dtauangle(j,1,1,i)=-1.0d0/sing*dcostau(j,1,1,i)
            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)
+           dtauangle(j,1,2,i)=-1.0d0/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)
+           dtauangle(j,1,3,i)=-1.0d0/sing*dcostau(j,1,3,i)
 c         write (iout,*) "else",i
          enddo
         endif
@@ -250,7 +255,7 @@ c        enddo
 c    Obtaining the gamma derivatives from sine derivative                                
        if (tauangle(2,i).gt.-pi4.and.tauangle(2,i).le.pi4.or.
      &     tauangle(2,i).gt.pi34.and.tauangle(2,i).le.pi.or.
-     &     tauangle(2,i).gt.-pi.and.tauangle(2,i).le.-pi34) then
+     &     tauangle(2,i).ge.-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_norm(1,i-1+nres),vp2)
          call vecpr(dc_norm(1,i-3),dc_norm(1,i-2),vp3)
@@ -282,15 +287,15 @@ c   Obtaining the gamma derivatives from cosine derivative
            dcostau(j,2,1,i)=fac1*dcostheta(j,1,i-1)+fac3*
      &     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)
+           dtauangle(j,2,1,i)=-1.0d0/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)
-           dtauangle(j,2,2,i)=-1/sing*dcostau(j,2,2,i)
+           dtauangle(j,2,2,i)=-1.0d0/sing*dcostau(j,2,2,i)
            dcostau(j,2,3,i)=fac2*dcosomicron(j,1,2,i)+fac4*
      &     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)
+           dtauangle(j,2,3,i)=-1.0d0/sing*dcostau(j,2,3,i)
 c        write(iout,*) i,j,"else", dtauangle(j,2,3,i) 
          enddo
         endif                                                                                            
@@ -313,10 +318,10 @@ c the conventional case
         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)
+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)
-        enddo
+C        enddo
         scalp=scalar(dc_norm2(1,i-2+nres),dc_norm(1,i-1+nres))
         fac0=1.0d0/(sint1*sint)
         fac1=cost*fac0
@@ -326,7 +331,7 @@ c        dc_norm2(j,i-1+nres)=-dc_norm(j,i-1+nres)
 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
+     &     tauangle(3,i).ge.-pi.and.tauangle(3,i).le.-pi34) then
          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)
@@ -355,15 +360,15 @@ c   Obtaining the gamma derivatives from cosine derivative
            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)
+           dtauangle(j,3,1,i)=-1.0d0/sing*dcostau(j,3,1,i)
            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)
+           dtauangle(j,3,2,i)=-1.0d0/sing*dcostau(j,3,2,i)
            dcostau(j,3,3,i)=fac2*dcosomicron(j,1,2,i)+fac4*
      &     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)
+           dtauangle(j,3,3,i)=-1.0d0/sing*dcostau(j,3,3,i)
 c          write(iout,*) "else",i 
          enddo
         endif                                                                                            
@@ -403,7 +408,7 @@ c   Derivatives of side-chain angles alpha and omega
 c obtaining the derivatives of omega from sines            
             if(omeg(i).gt.-pi4.and.omeg(i).le.pi4.or.
      &         omeg(i).gt.pi34.and.omeg(i).le.pi.or.
-     &         omeg(i).gt.-pi.and.omeg(i).le.-pi34) then
+     &         omeg(i).ge.-pi.and.omeg(i).le.-pi34) then
                fac15=dcos(theta(i+1))/(dsin(theta(i+1))*
      &        dsin(theta(i+1)))
                fac16=dcos(alph(i))/(dsin(alph(i))*dsin(alph(i)))