source/unres/src-HCD-5D/intcartderiv.F

   1       subroutine intcartderiv
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4 #ifdef MPI
   5       include 'mpif.h'
   6 #endif
   7       include 'COMMON.SETUP'
   8       include 'COMMON.CHAIN'
   9       include 'COMMON.VAR'
  10       include 'COMMON.GEO'
  11       include 'COMMON.INTERACT'
  12       include 'COMMON.DERIV'
  13       include 'COMMON.IOUNITS'
  14       include 'COMMON.LOCAL'
  15       include 'COMMON.SCCOR'
  16       double precision dcostheta(3,2,maxres),
  17      & dcosphi(3,3,maxres),dsinphi(3,3,maxres),
  18      & dcosalpha(3,3,maxres),dcosomega(3,3,maxres),
  19      & dsinomega(3,3,maxres),vo1(3),vo2(3),vo3(3),
  20      & dummy(3),vp1(3),vp2(3),vp3(3),vpp1(3),n(3)
  21
  22 #if defined(MPI) && defined(PARINTDER)
  23       if (nfgtasks.gt.1 .and. me.eq.king)
  24      &  call MPI_Bcast(8,1,MPI_INTEGER,king,FG_COMM,IERROR)
  25 #endif
  26       pi4 = 0.5d0*pipol
  27       pi34 = 3*pi4
  28
  29 c      write (iout,*) "iphi1_start",iphi1_start," iphi1_end",iphi1_end
  30       do i=1,nres
  31         do j=1,3
  32           dtheta(j,1,i)=0.0d0
  33           dtheta(j,2,i)=0.0d0
  34           dphi(j,1,i)=0.0d0
  35           dphi(j,2,i)=0.0d0
  36           dphi(j,3,i)=0.0d0
  37         enddo
  38       enddo
  39 c Derivatives of thetas
  40 #if defined(MPI) && defined(PARINTDER)
  41 c We need dtheta(:,:,i-1) to compute dphi(:,:,i)
  42       do i=max0(ithet_start-1,3),ithet_end
  43 #else
  44       do i=3,nres
  45 #endif
  46         cost=dcos(theta(i))
  47         sint=sqrt(1-cost*cost)
  48         do j=1,3
  49           dcostheta(j,1,i)=-(dc_norm(j,i-1)+cost*dc_norm(j,i-2))/
  50      &    vbld(i-1)
  51 c          if (itype(i-1).ne.ntyp1)
  52           dtheta(j,1,i)=-dcostheta(j,1,i)/sint
  53           dcostheta(j,2,i)=-(dc_norm(j,i-2)+cost*dc_norm(j,i-1))/
  54      &    vbld(i)
  55 c          if (itype(i-1).ne.ntyp1)
  56           dtheta(j,2,i)=-dcostheta(j,2,i)/sint
  57         enddo
  58       enddo
  59 #if defined(MPI) && defined(PARINTDER)
  60 c We need dtheta(:,:,i-1) to compute dphi(:,:,i)
  61       do i=max0(ithet_start-1,3),ithet_end
  62 #else
  63       do i=3,nres
  64 #endif
  65       if ((itype(i-1).ne.10).and.(itype(i-1).ne.ntyp1)) then
  66         cost1=dcos(omicron(1,i))
  67         sint1=sqrt(1-cost1*cost1)
  68         cost2=dcos(omicron(2,i))
  69         sint2=sqrt(1-cost2*cost2)
  70        do j=1,3
  71 CC Calculate derivative over first omicron (Cai-2,Cai-1,SCi-1)
  72           dcosomicron(j,1,1,i)=-(dc_norm(j,i-1+nres)+
  73      &    cost1*dc_norm(j,i-2))/
  74      &    vbld(i-1)
  75           domicron(j,1,1,i)=-1/sint1*dcosomicron(j,1,1,i)
  76           dcosomicron(j,1,2,i)=-(dc_norm(j,i-2)
  77      &    +cost1*(dc_norm(j,i-1+nres)))/
  78      &    vbld(i-1+nres)
  79           domicron(j,1,2,i)=-1/sint1*dcosomicron(j,1,2,i)
  80 CC Calculate derivative over second omicron Sci-1,Cai-1 Cai
  81 CC Looks messy but better than if in loop
  82           dcosomicron(j,2,1,i)=-(-dc_norm(j,i-1+nres)
  83      &    +cost2*dc_norm(j,i-1))/
  84      &    vbld(i)
  85           domicron(j,2,1,i)=-1/sint2*dcosomicron(j,2,1,i)
  86           dcosomicron(j,2,2,i)=-(dc_norm(j,i-1)
  87      &     +cost2*(-dc_norm(j,i-1+nres)))/
  88      &    vbld(i-1+nres)
  89 c          write(iout,*) "vbld", i,itype(i),vbld(i-1+nres)
  90           domicron(j,2,2,i)=-1/sint2*dcosomicron(j,2,2,i)
  91         enddo
  92        endif
  93       enddo
  94
  95 c Derivatives of phi:
  96 c If phi is 0 or 180 degrees, then the formulas
  97 c have to be derived by power series expansion of the
  98 c conventional formulas around 0 and 180.
  99 #ifdef PARINTDER
 100       do i=iphi1_start,iphi1_end
 101 #else
 102       do i=4,nres
 103 #endif
 104 c        if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1
 105 c     &      .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle
 106 c the conventional case
 107         sint=dsin(theta(i))
 108         sint1=dsin(theta(i-1))
 109         sing=dsin(phi(i))
 110         cost=dcos(theta(i))
 111         cost1=dcos(theta(i-1))
 112         cosg=dcos(phi(i))
 113         scalp=scalar(dc_norm(1,i-3),dc_norm(1,i-1))
 114         fac0=1.0d0/(sint1*sint)
 115         fac1=cost*fac0
 116         fac2=cost1*fac0
 117         fac3=cosg*cost1/(sint1*sint1)
 118         fac4=cosg*cost/(sint*sint)
 119 c    Obtaining the gamma derivatives from sine derivative
 120        if (phi(i).gt.-pi4.and.phi(i).le.pi4.or.
 121      &     phi(i).gt.pi34.and.phi(i).le.pi.or.
 122      &     phi(i).ge.-pi.and.phi(i).le.-pi34) then
 123          call vecpr(dc_norm(1,i-1),dc_norm(1,i-2),vp1)
 124          call vecpr(dc_norm(1,i-3),dc_norm(1,i-1),vp2)
 125          call vecpr(dc_norm(1,i-3),dc_norm(1,i-2),vp3)
 126          do j=1,3
 127             ctgt=cost/sint
 128             ctgt1=cost1/sint1
 129             cosg_inv=1.0d0/cosg
 130 c            if (itype(i-1).ne.ntyp1 .and. itype(i-2).ne.ntyp1) then
 131       dsinphi(j,1,i)=-sing*ctgt1*dtheta(j,1,i-1)
 132      &        -(fac0*vp1(j)+sing*dc_norm(j,i-3))*vbld_inv(i-2)
 133             dphi(j,1,i)=cosg_inv*dsinphi(j,1,i)
 134             dsinphi(j,2,i)=
 135      &        -sing*(ctgt1*dtheta(j,2,i-1)+ctgt*dtheta(j,1,i))
 136      &        -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
 137             dphi(j,2,i)=cosg_inv*dsinphi(j,2,i)
 138             dsinphi(j,3,i)=-sing*ctgt*dtheta(j,2,i)
 139      &        +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i)
 140 c     &        +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
 141             dphi(j,3,i)=cosg_inv*dsinphi(j,3,i)
 142 c            endif
 143 c Bug fixed 3/24/05 (AL)
 144          enddo
 145 c   Obtaining the gamma derivatives from cosine derivative
 146         else
 147            do j=1,3
 148 c           if (itype(i-1).ne.ntyp1 .and. itype(i-2).ne.ntyp1) then
 149            dcosphi(j,1,i)=fac1*dcostheta(j,1,i-1)+fac3*
 150      &     dcostheta(j,1,i-1)-fac0*(dc_norm(j,i-1)-scalp*
 151      &     dc_norm(j,i-3))/vbld(i-2)
 152            dphi(j,1,i)=-1/sing*dcosphi(j,1,i)
 153            dcosphi(j,2,i)=fac1*dcostheta(j,2,i-1)+fac2*
 154      &     dcostheta(j,1,i)+fac3*dcostheta(j,2,i-1)+fac4*
 155      &     dcostheta(j,1,i)
 156            dphi(j,2,i)=-1/sing*dcosphi(j,2,i)
 157            dcosphi(j,3,i)=fac2*dcostheta(j,2,i)+fac4*
 158      &     dcostheta(j,2,i)-fac0*(dc_norm(j,i-3)-scalp*
 159      &     dc_norm(j,i-1))/vbld(i)
 160            dphi(j,3,i)=-1/sing*dcosphi(j,3,i)
 161 c           endif
 162          enddo
 163         endif
 164       enddo
 165 Calculate derivative of Tauangle
 166       do i=1,nres-1
 167        do j=1,3
 168         dc_norm2(j,i+nres)=-dc_norm(j,i+nres)
 169        enddo
 170       enddo
 171 #ifdef PARINTDER
 172       do i=itau_start,itau_end
 173 #else
 174       do i=3,nres
 175 #endif
 176        if ((itype(i-2).eq.ntyp1).or.(itype(i-2).eq.10)) cycle
 177 c       if ((itype(i-2).eq.ntyp1).or.(itype(i-2).eq.10).or.
 178 c     &     (itype(i-1).eq.ntyp1).or.(itype(i).eq.ntyp1)) cycle
 179 cc dtauangle(j,intertyp,dervityp,residue number)
 180 cc INTERTYP=1 SC...Ca...Ca..Ca
 181 c the conventional case
 182         sint=dsin(theta(i))
 183         sint1=dsin(omicron(2,i-1))
 184         sing=dsin(tauangle(1,i))
 185         cost=dcos(theta(i))
 186         cost1=dcos(omicron(2,i-1))
 187         cosg=dcos(tauangle(1,i))
 188         do j=1,3
 189         dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres)
 190 cc       write(iout,*) dc_norm2(j,i-2+nres),"dcnorm"
 191         enddo
 192         scalp=scalar(dc_norm2(1,i-2+nres),dc_norm(1,i-1))
 193         fac0=1.0d0/(sint1*sint)
 194         fac1=cost*fac0
 195         fac2=cost1*fac0
 196         fac3=cosg*cost1/(sint1*sint1)
 197         fac4=cosg*cost/(sint*sint)
 198 cc         write(iout,*) "faki",fac0,fac1,fac2,fac3,fac4
 199 c    Obtaining the gamma derivatives from sine derivative
 200        if (tauangle(1,i).gt.-pi4.and.tauangle(1,i).le.pi4.or.
 201      &     tauangle(1,i).gt.pi34.and.tauangle(1,i).le.pi.or.
 202      &     tauangle(1,i).ge.-pi.and.tauangle(1,i).le.-pi34) then
 203          call vecpr(dc_norm(1,i-1),dc_norm(1,i-2),vp1)
 204          call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-1),vp2)
 205          call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-2),vp3)
 206         do j=1,3
 207             ctgt=cost/sint
 208             ctgt1=cost1/sint1
 209             cosg_inv=1.0d0/cosg
 210             dsintau(j,1,1,i)=-sing*ctgt1*domicron(j,2,2,i-1)
 211      &-(fac0*vp1(j)+sing*(dc_norm2(j,i-2+nres)))
 212      & *vbld_inv(i-2+nres)
 213             dtauangle(j,1,1,i)=cosg_inv*dsintau(j,1,1,i)
 214             dsintau(j,1,2,i)=
 215      &        -sing*(ctgt1*domicron(j,2,1,i-1)+ctgt*dtheta(j,1,i))
 216      &        -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
 217 c            write(iout,*) "dsintau", dsintau(j,1,2,i)
 218             dtauangle(j,1,2,i)=cosg_inv*dsintau(j,1,2,i)
 219 c Bug fixed 3/24/05 (AL)
 220             dsintau(j,1,3,i)=-sing*ctgt*dtheta(j,2,i)
 221      &        +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i)
 222 c     &        +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
 223             dtauangle(j,1,3,i)=cosg_inv*dsintau(j,1,3,i)
 224          enddo
 225 c   Obtaining the gamma derivatives from cosine derivative
 226         else
 227            do j=1,3
 228            dcostau(j,1,1,i)=fac1*dcosomicron(j,2,2,i-1)+fac3*
 229      &     dcosomicron(j,2,2,i-1)-fac0*(dc_norm(j,i-1)-scalp*
 230      &     (dc_norm2(j,i-2+nres)))/vbld(i-2+nres)
 231            dtauangle(j,1,1,i)=-1/sing*dcostau(j,1,1,i)
 232            dcostau(j,1,2,i)=fac1*dcosomicron(j,2,1,i-1)+fac2*
 233      &     dcostheta(j,1,i)+fac3*dcosomicron(j,2,1,i-1)+fac4*
 234      &     dcostheta(j,1,i)
 235            dtauangle(j,1,2,i)=-1/sing*dcostau(j,1,2,i)
 236            dcostau(j,1,3,i)=fac2*dcostheta(j,2,i)+fac4*
 237      &     dcostheta(j,2,i)-fac0*(-dc_norm(j,i-2+nres)-scalp*
 238      &     dc_norm(j,i-1))/vbld(i)
 239            dtauangle(j,1,3,i)=-1/sing*dcostau(j,1,3,i)
 240 c         write (iout,*) "else",i
 241          enddo
 242         endif
 243 c        do k=1,3
 244 c        write(iout,*) "tu",i,k,(dtauangle(j,1,k,i),j=1,3)
 245 c        enddo
 246       enddo
 247 CC Second case Ca...Ca...Ca...SC
 248 #ifdef PARINTDER
 249       do i=itau_start,itau_end
 250 #else
 251       do i=4,nres
 252 #endif
 253        if ((itype(i-1).eq.ntyp1).or.(itype(i-1).eq.10).or.
 254      &    (itype(i-2).eq.ntyp1).or.(itype(i-3).eq.ntyp1)) cycle
 255 c the conventional case
 256         sint=dsin(omicron(1,i))
 257         sint1=dsin(theta(i-1))
 258         sing=dsin(tauangle(2,i))
 259         cost=dcos(omicron(1,i))
 260         cost1=dcos(theta(i-1))
 261         cosg=dcos(tauangle(2,i))
 262 c        do j=1,3
 263 c        dc_norm2(j,i-1+nres)=-dc_norm(j,i-1+nres)
 264 c        enddo
 265         scalp=scalar(dc_norm(1,i-3),dc_norm(1,i-1+nres))
 266         fac0=1.0d0/(sint1*sint)
 267         fac1=cost*fac0
 268         fac2=cost1*fac0
 269         fac3=cosg*cost1/(sint1*sint1)
 270         fac4=cosg*cost/(sint*sint)
 271 c    Obtaining the gamma derivatives from sine derivative
 272        if (tauangle(2,i).gt.-pi4.and.tauangle(2,i).le.pi4.or.
 273      &     tauangle(2,i).gt.pi34.and.tauangle(2,i).le.pi.or.
 274      &     tauangle(2,i).gt.-pi.and.tauangle(2,i).le.-pi34) then
 275          call vecpr(dc_norm2(1,i-1+nres),dc_norm(1,i-2),vp1)
 276          call vecpr(dc_norm(1,i-3),dc_norm(1,i-1+nres),vp2)
 277          call vecpr(dc_norm(1,i-3),dc_norm(1,i-2),vp3)
 278         do j=1,3
 279             ctgt=cost/sint
 280             ctgt1=cost1/sint1
 281             cosg_inv=1.0d0/cosg
 282             dsintau(j,2,1,i)=-sing*ctgt1*dtheta(j,1,i-1)
 283      &        +(fac0*vp1(j)-sing*dc_norm(j,i-3))*vbld_inv(i-2)
 284 c       write(iout,*) i,j,dsintau(j,2,1,i),sing*ctgt1*dtheta(j,1,i-1),
 285 c     &fac0*vp1(j),sing*dc_norm(j,i-3),vbld_inv(i-2),"dsintau(2,1)"
 286             dtauangle(j,2,1,i)=cosg_inv*dsintau(j,2,1,i)
 287             dsintau(j,2,2,i)=
 288      &        -sing*(ctgt1*dtheta(j,2,i-1)+ctgt*domicron(j,1,1,i))
 289      &        -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
 290 c            write(iout,*) "sprawdzenie",i,j,sing*ctgt1*dtheta(j,2,i-1),
 291 c     & sing*ctgt*domicron(j,1,2,i),
 292 c     & (fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
 293             dtauangle(j,2,2,i)=cosg_inv*dsintau(j,2,2,i)
 294 c Bug fixed 3/24/05 (AL)
 295             dsintau(j,2,3,i)=-sing*ctgt*domicron(j,1,2,i)
 296      &       +(fac0*vp3(j)-sing*dc_norm(j,i-1+nres))*vbld_inv(i-1+nres)
 297 c     &        +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
 298             dtauangle(j,2,3,i)=cosg_inv*dsintau(j,2,3,i)
 299          enddo
 300 c   Obtaining the gamma derivatives from cosine derivative
 301         else
 302            do j=1,3
 303            dcostau(j,2,1,i)=fac1*dcostheta(j,1,i-1)+fac3*
 304      &     dcostheta(j,1,i-1)-fac0*(dc_norm(j,i-1+nres)-scalp*
 305      &     dc_norm(j,i-3))/vbld(i-2)
 306            dtauangle(j,2,1,i)=-1/sing*dcostau(j,2,1,i)
 307            dcostau(j,2,2,i)=fac1*dcostheta(j,2,i-1)+fac2*
 308      &     dcosomicron(j,1,1,i)+fac3*dcostheta(j,2,i-1)+fac4*
 309      &     dcosomicron(j,1,1,i)
 310            dtauangle(j,2,2,i)=-1/sing*dcostau(j,2,2,i)
 311            dcostau(j,2,3,i)=fac2*dcosomicron(j,1,2,i)+fac4*
 312      &     dcosomicron(j,1,2,i)-fac0*(dc_norm(j,i-3)-scalp*
 313      &     dc_norm(j,i-1+nres))/vbld(i-1+nres)
 314            dtauangle(j,2,3,i)=-1/sing*dcostau(j,2,3,i)
 315 c        write(iout,*) i,j,"else", dtauangle(j,2,3,i)
 316          enddo
 317         endif
 318       enddo
 319
 320 CCC third case SC...Ca...Ca...SC
 321 #ifdef PARINTDER
 322
 323       do i=itau_start,itau_end
 324 #else
 325       do i=3,nres
 326 #endif
 327 c the conventional case
 328       if ((itype(i-1).eq.ntyp1).or.(itype(i-1).eq.10).or.
 329      &(itype(i-2).eq.ntyp1).or.(itype(i-2).eq.10)) cycle
 330         sint=dsin(omicron(1,i))
 331         sint1=dsin(omicron(2,i-1))
 332         sing=dsin(tauangle(3,i))
 333         cost=dcos(omicron(1,i))
 334         cost1=dcos(omicron(2,i-1))
 335         cosg=dcos(tauangle(3,i))
 336 c        write (iout,*) "i",i," omicron",omicron(1,i),omicron(2,i)
 337 c        write (iout,*) "i",i," tauangle",tauangle(1,i),tauangle(2,i),
 338 c     &    tauangle(3,i)
 339         do j=1,3
 340         dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres)
 341 c        dc_norm2(j,i-1+nres)=-dc_norm(j,i-1+nres)
 342         enddo
 343         scalp=scalar(dc_norm2(1,i-2+nres),dc_norm(1,i-1+nres))
 344         fac0=1.0d0/(sint1*sint)
 345         fac1=cost*fac0
 346         fac2=cost1*fac0
 347         fac3=cosg*cost1/(sint1*sint1)
 348         fac4=cosg*cost/(sint*sint)
 349 c    Obtaining the gamma derivatives from sine derivative
 350        if (tauangle(3,i).gt.-pi4.and.tauangle(3,i).le.pi4.or.
 351      &     tauangle(3,i).gt.pi34.and.tauangle(3,i).le.pi.or.
 352      &     tauangle(3,i).ge.-pi.and.tauangle(3,i).le.-pi34) then
 353          call vecpr(dc_norm(1,i-1+nres),dc_norm(1,i-2),vp1)
 354          call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-1+nres),vp2)
 355          call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-2),vp3)
 356         do j=1,3
 357             ctgt=cost/sint
 358             ctgt1=cost1/sint1
 359             cosg_inv=1.0d0/cosg
 360             dsintau(j,3,1,i)=-sing*ctgt1*domicron(j,2,2,i-1)
 361      &        -(fac0*vp1(j)-sing*dc_norm(j,i-2+nres))
 362      &        *vbld_inv(i-2+nres)
 363             dtauangle(j,3,1,i)=cosg_inv*dsintau(j,3,1,i)
 364             dsintau(j,3,2,i)=
 365      &        -sing*(ctgt1*domicron(j,2,1,i-1)+ctgt*domicron(j,1,1,i))
 366      &        -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
 367             dtauangle(j,3,2,i)=cosg_inv*dsintau(j,3,2,i)
 368 c Bug fixed 3/24/05 (AL)
 369             dsintau(j,3,3,i)=-sing*ctgt*domicron(j,1,2,i)
 370      &        +(fac0*vp3(j)-sing*dc_norm(j,i-1+nres))
 371      &        *vbld_inv(i-1+nres)
 372 c     &        +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
 373             dtauangle(j,3,3,i)=cosg_inv*dsintau(j,3,3,i)
 374          enddo
 375 c   Obtaining the gamma derivatives from cosine derivative
 376         else
 377            do j=1,3
 378            dcostau(j,3,1,i)=fac1*dcosomicron(j,2,2,i-1)+fac3*
 379      &     dcosomicron(j,2,2,i-1)-fac0*(dc_norm(j,i-1+nres)-scalp*
 380      &     dc_norm2(j,i-2+nres))/vbld(i-2+nres)
 381            dtauangle(j,3,1,i)=-1/sing*dcostau(j,3,1,i)
 382            dcostau(j,3,2,i)=fac1*dcosomicron(j,2,1,i-1)+fac2*
 383      &     dcosomicron(j,1,1,i)+fac3*dcosomicron(j,2,1,i-1)+fac4*
 384      &     dcosomicron(j,1,1,i)
 385            dtauangle(j,3,2,i)=-1/sing*dcostau(j,3,2,i)
 386            dcostau(j,3,3,i)=fac2*dcosomicron(j,1,2,i)+fac4*
 387      &     dcosomicron(j,1,2,i)-fac0*(dc_norm2(j,i-2+nres)-scalp*
 388      &     dc_norm(j,i-1+nres))/vbld(i-1+nres)
 389            dtauangle(j,3,3,i)=-1/sing*dcostau(j,3,3,i)
 390 c          write(iout,*) "else",i
 391          enddo
 392         endif
 393       enddo
 394
 395 #ifdef CRYST_SC
 396 c   Derivatives of side-chain angles alpha and omega
 397 #if defined(MPI) && defined(PARINTDER)
 398         do i=ibond_start,ibond_end
 399 #else
 400         do i=2,nres-1
 401 #endif
 402           if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 403              fac5=1.0d0/dsqrt(2*(1+dcos(theta(i+1))))
 404              fac6=fac5/vbld(i)
 405              fac7=fac5*fac5
 406              fac8=fac5/vbld(i+1)
 407              fac9=fac5/vbld(i+nres)
 408              scala1=scalar(dc_norm(1,i-1),dc_norm(1,i+nres))
 409              scala2=scalar(dc_norm(1,i),dc_norm(1,i+nres))
 410              cosa=dsqrt(0.5d0/(1.0d0+dcos(theta(i+1))))*(
 411      &       scalar(dC_norm(1,i),dC_norm(1,i+nres))
 412      &       -scalar(dC_norm(1,i-1),dC_norm(1,i+nres)))
 413              sina=sqrt(1-cosa*cosa)
 414              sino=dsin(omeg(i))
 415 c             write (iout,*) "i",i," cosa",cosa," sina",sina," sino",sino
 416              do j=1,3
 417                 dcosalpha(j,1,i)=fac6*(scala1*dc_norm(j,i-1)-
 418      &          dc_norm(j,i+nres))-cosa*fac7*dcostheta(j,1,i+1)
 419                 dalpha(j,1,i)=-1/sina*dcosalpha(j,1,i)
 420                 dcosalpha(j,2,i)=fac8*(dc_norm(j,i+nres)-
 421      &          scala2*dc_norm(j,i))-cosa*fac7*dcostheta(j,2,i+1)
 422                 dalpha(j,2,i)=-1/sina*dcosalpha(j,2,i)
 423                 dcosalpha(j,3,i)=(fac9*(dc_norm(j,i)-
 424      &          dc_norm(j,i-1))-(cosa*dc_norm(j,i+nres))/
 425      &          vbld(i+nres))
 426                 dalpha(j,3,i)=-1/sina*dcosalpha(j,3,i)
 427             enddo
 428 c obtaining the derivatives of omega from sines
 429             if(omeg(i).gt.-pi4.and.omeg(i).le.pi4.or.
 430      &         omeg(i).gt.pi34.and.omeg(i).le.pi.or.
 431      &         omeg(i).ge.-pi.and.omeg(i).le.-pi34) then
 432                fac15=dcos(theta(i+1))/(dsin(theta(i+1))*
 433      &         dsin(theta(i+1)))
 434                fac16=dcos(alph(i))/(dsin(alph(i))*dsin(alph(i)))
 435                fac17=1.0d0/(dsin(theta(i+1))*dsin(alph(i)))
 436                call vecpr(dc_norm(1,i+nres),dc_norm(1,i),vo1)
 437                call vecpr(dc_norm(1,i+nres),dc_norm(1,i-1),vo2)
 438                call vecpr(dc_norm(1,i),dc_norm(1,i-1),vo3)
 439                coso_inv=1.0d0/dcos(omeg(i))
 440                do j=1,3
 441                  dsinomega(j,1,i)=sino*(fac15*dcostheta(j,1,i+1)
 442      &           +fac16*dcosalpha(j,1,i))-fac17/vbld(i)*vo1(j)-(
 443      &           sino*dc_norm(j,i-1))/vbld(i)
 444                  domega(j,1,i)=coso_inv*dsinomega(j,1,i)
 445                  dsinomega(j,2,i)=sino*(fac15*dcostheta(j,2,i+1)
 446      &           +fac16*dcosalpha(j,2,i))+fac17/vbld(i+1)*vo2(j)
 447      &           -sino*dc_norm(j,i)/vbld(i+1)
 448                  domega(j,2,i)=coso_inv*dsinomega(j,2,i)
 449                  dsinomega(j,3,i)=sino*fac16*dcosalpha(j,3,i)-
 450      &           fac17/vbld(i+nres)*vo3(j)-sino*dc_norm(j,i+nres)/
 451      &           vbld(i+nres)
 452                  domega(j,3,i)=coso_inv*dsinomega(j,3,i)
 453               enddo
 454            else
 455 c   obtaining the derivatives of omega from cosines
 456              fac10=sqrt(0.5d0*(1-dcos(theta(i+1))))
 457              fac11=sqrt(0.5d0*(1+dcos(theta(i+1))))
 458              fac12=fac10*sina
 459              fac13=fac12*fac12
 460              fac14=sina*sina
 461              do j=1,3
 462                 dcosomega(j,1,i)=(-(0.25d0*cosa/fac11*
 463      &          dcostheta(j,1,i+1)+fac11*dcosalpha(j,1,i))*fac12+
 464      &          (0.25d0/fac10*sina*dcostheta(j,1,i+1)+cosa/sina*
 465      &          fac10*dcosalpha(j,1,i))*(scala2-fac11*cosa))/fac13
 466                 domega(j,1,i)=-1/sino*dcosomega(j,1,i)
 467                 dcosomega(j,2,i)=(((dc_norm(j,i+nres)-scala2*
 468      &          dc_norm(j,i))/vbld(i+1)-0.25d0*cosa/fac11*
 469      &          dcostheta(j,2,i+1)-fac11*dcosalpha(j,2,i))*fac12+
 470      &          (scala2-fac11*cosa)*(0.25d0*sina/fac10*
 471      &          dcostheta(j,2,i+1)+fac10*cosa/sina*dcosalpha(j,2,i)
 472      &          ))/fac13
 473                 domega(j,2,i)=-1/sino*dcosomega(j,2,i)
 474                 dcosomega(j,3,i)=1/fac10*((1/vbld(i+nres)*(dc_norm(j,i)-
 475      &          scala2*dc_norm(j,i+nres))-fac11*dcosalpha(j,3,i))*sina+
 476      &          (scala2-fac11*cosa)*(cosa/sina*dcosalpha(j,3,i)))/fac14
 477                 domega(j,3,i)=-1/sino*dcosomega(j,3,i)
 478             enddo
 479           endif
 480          else
 481            do j=1,3
 482              do k=1,3
 483                dalpha(k,j,i)=0.0d0
 484                domega(k,j,i)=0.0d0
 485              enddo
 486            enddo
 487          endif
 488        enddo
 489 #endif
 490 #if defined(MPI) && defined(PARINTDER)
 491       if (nfgtasks.gt.1) then
 492 #ifdef DEBUG
 493 cd      write (iout,*) "Gather dtheta"
 494 cd      call flush(iout)
 495       write (iout,*) "dtheta before gather"
 496       do i=1,nres
 497         write (iout,'(i3,3(3f8.5,3x))') i,((dtheta(j,k,i),k=1,3),j=1,2)
 498       enddo
 499 #endif
 500       call MPI_Gatherv(dtheta(1,1,ithet_start),ithet_count(fg_rank),
 501      &  MPI_THET,dtheta(1,1,1),ithet_count(0),ithet_displ(0),MPI_THET,
 502      &  king,FG_COMM,IERROR)
 503 #ifdef DEBUG
 504 cd      write (iout,*) "Gather dphi"
 505 cd      call flush(iout)
 506       write (iout,*) "dphi before gather"
 507       do i=1,nres
 508         write (iout,'(i3,3(3f8.5,3x))') i,((dphi(j,k,i),k=1,3),j=1,3)
 509       enddo
 510 #endif
 511       call MPI_Gatherv(dphi(1,1,iphi1_start),iphi1_count(fg_rank),
 512      &  MPI_GAM,dphi(1,1,1),iphi1_count(0),iphi1_displ(0),MPI_GAM,
 513      &  king,FG_COMM,IERROR)
 514 cd      write (iout,*) "Gather dalpha"
 515 cd      call flush(iout)
 516 #ifdef CRYST_SC
 517       call MPI_Gatherv(dalpha(1,1,ibond_start),ibond_count(fg_rank),
 518      &  MPI_GAM,dalpha(1,1,1),ibond_count(0),ibond_displ(0),MPI_GAM,
 519      &  king,FG_COMM,IERROR)
 520 cd      write (iout,*) "Gather domega"
 521 cd      call flush(iout)
 522       call MPI_Gatherv(domega(1,1,ibond_start),ibond_count(fg_rank),
 523      &  MPI_GAM,domega(1,1,1),ibond_count(0),ibond_displ(0),MPI_GAM,
 524      &  king,FG_COMM,IERROR)
 525 #endif
 526       endif
 527 #endif
 528 #ifdef DEBUG
 529       write (iout,*) "dtheta after gather"
 530       do i=1,nres
 531         write (iout,'(i3,3(3f8.5,3x))') i,((dtheta(j,k,i),j=1,3),k=1,2)
 532       enddo
 533       write (iout,*) "dphi after gather"
 534       do i=1,nres
 535         write (iout,'(i3,3(3f8.5,3x))') i,((dphi(j,k,i),j=1,3),k=1,3)
 536       enddo
 537       write (iout,*) "dalpha after gather"
 538       do i=1,nres
 539         write (iout,'(i3,3(3f8.5,3x))') i,((dalpha(j,k,i),j=1,3),k=1,3)
 540       enddo
 541       write (iout,*) "domega after gather"
 542       do i=1,nres
 543         write (iout,'(i3,3(3f8.5,3x))') i,((domega(j,k,i),j=1,3),k=1,3)
 544       enddo
 545       write (iout,*) "dtauangle after gather"
 546       do i=1,nres
 547         write (iout,'(i3,3(3f8.5,3x))') i,
 548      &         (((dtauangle(j,k,l,i),j=1,3),k=1,3),l=1,3)
 549       enddo
 550 #endif
 551       return
 552       end
 553
 554       subroutine checkintcartgrad
 555       implicit real*8 (a-h,o-z)
 556       include 'DIMENSIONS'
 557 #ifdef MPI
 558       include 'mpif.h'
 559 #endif
 560       include 'COMMON.CHAIN'
 561       include 'COMMON.VAR'
 562       include 'COMMON.GEO'
 563       include 'COMMON.INTERACT'
 564       include 'COMMON.DERIV'
 565       include 'COMMON.IOUNITS'
 566       include 'COMMON.SETUP'
 567       double precision dthetanum(3,2,maxres),dphinum(3,3,maxres)
 568      & ,dalphanum(3,3,maxres), domeganum(3,3,maxres)
 569       double precision theta_s(maxres),phi_s(maxres),alph_s(maxres),
 570      & omeg_s(maxres),dc_norm_s(3)
 571       double precision aincr /1.0d-5/
 572
 573       do i=1,nres
 574         phi_s(i)=phi(i)
 575         theta_s(i)=theta(i)
 576         alph_s(i)=alph(i)
 577         omeg_s(i)=omeg(i)
 578       enddo
 579 c Check theta gradient
 580       write (iout,*)
 581      & "Analytical (upper) and numerical (lower) gradient of theta"
 582       write (iout,*)
 583       do i=3,nres
 584         do j=1,3
 585           dcji=dc(j,i-2)
 586           dc(j,i-2)=dcji+aincr
 587           call chainbuild_cart
 588           call int_from_cart1(.false.)
 589           dthetanum(j,1,i)=(theta(i)-theta_s(i))/aincr
 590           dc(j,i-2)=dcji
 591           dcji=dc(j,i-1)
 592           dc(j,i-1)=dc(j,i-1)+aincr
 593           call chainbuild_cart
 594           dthetanum(j,2,i)=(theta(i)-theta_s(i))/aincr
 595           dc(j,i-1)=dcji
 596         enddo
 597         write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dtheta(j,1,i),j=1,3),
 598      &    (dtheta(j,2,i),j=1,3)
 599         write (iout,'(5x,3f10.5,5x,3f10.5)') (dthetanum(j,1,i),j=1,3),
 600      &    (dthetanum(j,2,i),j=1,3)
 601         write (iout,'(5x,3f10.5,5x,3f10.5)')
 602      &    (dthetanum(j,1,i)/dtheta(j,1,i),j=1,3),
 603      &    (dthetanum(j,2,i)/dtheta(j,2,i),j=1,3)
 604         write (iout,*)
 605       enddo
 606 c Check gamma gradient
 607       write (iout,*)
 608      & "Analytical (upper) and numerical (lower) gradient of gamma"
 609       do i=4,nres
 610         do j=1,3
 611           dcji=dc(j,i-3)
 612           dc(j,i-3)=dcji+aincr
 613           call chainbuild_cart
 614           dphinum(j,1,i)=(phi(i)-phi_s(i))/aincr
 615           dc(j,i-3)=dcji
 616           dcji=dc(j,i-2)
 617           dc(j,i-2)=dcji+aincr
 618           call chainbuild_cart
 619           dphinum(j,2,i)=(phi(i)-phi_s(i))/aincr
 620           dc(j,i-2)=dcji
 621           dcji=dc(j,i-1)
 622           dc(j,i-1)=dc(j,i-1)+aincr
 623           call chainbuild_cart
 624           dphinum(j,3,i)=(phi(i)-phi_s(i))/aincr
 625           dc(j,i-1)=dcji
 626         enddo
 627         write (iout,'(i5,3(3f10.5,5x))') i,(dphi(j,1,i),j=1,3),
 628      &    (dphi(j,2,i),j=1,3),(dphi(j,3,i),j=1,3)
 629         write (iout,'(5x,3(3f10.5,5x))') (dphinum(j,1,i),j=1,3),
 630      &    (dphinum(j,2,i),j=1,3),(dphinum(j,3,i),j=1,3)
 631         write (iout,'(5x,3(3f10.5,5x))')
 632      &    (dphinum(j,1,i)/dphi(j,1,i),j=1,3),
 633      &    (dphinum(j,2,i)/dphi(j,2,i),j=1,3),
 634      &    (dphinum(j,3,i)/dphi(j,3,i),j=1,3)
 635         write (iout,*)
 636       enddo
 637 c Check alpha gradient
 638       write (iout,*)
 639      & "Analytical (upper) and numerical (lower) gradient of alpha"
 640       do i=2,nres-1
 641        if(itype(i).ne.10) then
 642             do j=1,3
 643               dcji=dc(j,i-1)
 644               dc(j,i-1)=dcji+aincr
 645               call chainbuild_cart
 646               dalphanum(j,1,i)=(alph(i)-alph_s(i))
 647      &        /aincr
 648               dc(j,i-1)=dcji
 649               dcji=dc(j,i)
 650               dc(j,i)=dcji+aincr
 651               call chainbuild_cart
 652               dalphanum(j,2,i)=(alph(i)-alph_s(i))
 653      &        /aincr
 654               dc(j,i)=dcji
 655               dcji=dc(j,i+nres)
 656               dc(j,i+nres)=dc(j,i+nres)+aincr
 657               call chainbuild_cart
 658               dalphanum(j,3,i)=(alph(i)-alph_s(i))
 659      &        /aincr
 660              dc(j,i+nres)=dcji
 661             enddo
 662           endif
 663         write (iout,'(i5,3(3f10.5,5x))') i,(dalpha(j,1,i),j=1,3),
 664      &    (dalpha(j,2,i),j=1,3),(dalpha(j,3,i),j=1,3)
 665         write (iout,'(5x,3(3f10.5,5x))') (dalphanum(j,1,i),j=1,3),
 666      &    (dalphanum(j,2,i),j=1,3),(dalphanum(j,3,i),j=1,3)
 667         write (iout,'(5x,3(3f10.5,5x))')
 668      &    (dalphanum(j,1,i)/dalpha(j,1,i),j=1,3),
 669      &    (dalphanum(j,2,i)/dalpha(j,2,i),j=1,3),
 670      &    (dalphanum(j,3,i)/dalpha(j,3,i),j=1,3)
 671         write (iout,*)
 672       enddo
 673 c     Check omega gradient
 674       write (iout,*)
 675      & "Analytical (upper) and numerical (lower) gradient of omega"
 676       do i=2,nres-1
 677        if(itype(i).ne.10) then
 678             do j=1,3
 679               dcji=dc(j,i-1)
 680               dc(j,i-1)=dcji+aincr
 681               call chainbuild_cart
 682               domeganum(j,1,i)=(omeg(i)-omeg_s(i))
 683      &        /aincr
 684               dc(j,i-1)=dcji
 685               dcji=dc(j,i)
 686               dc(j,i)=dcji+aincr
 687               call chainbuild_cart
 688               domeganum(j,2,i)=(omeg(i)-omeg_s(i))
 689      &        /aincr
 690               dc(j,i)=dcji
 691               dcji=dc(j,i+nres)
 692               dc(j,i+nres)=dc(j,i+nres)+aincr
 693               call chainbuild_cart
 694               domeganum(j,3,i)=(omeg(i)-omeg_s(i))
 695      &        /aincr
 696              dc(j,i+nres)=dcji
 697             enddo
 698           endif
 699         write (iout,'(i5,3(3f10.5,5x))') i,(domega(j,1,i),j=1,3),
 700      &    (domega(j,2,i),j=1,3),(domega(j,3,i),j=1,3)
 701         write (iout,'(5x,3(3f10.5,5x))') (domeganum(j,1,i),j=1,3),
 702      &    (domeganum(j,2,i),j=1,3),(domeganum(j,3,i),j=1,3)
 703         write (iout,'(5x,3(3f10.5,5x))')
 704      &    (domeganum(j,1,i)/domega(j,1,i),j=1,3),
 705      &    (domeganum(j,2,i)/domega(j,2,i),j=1,3),
 706      &    (domeganum(j,3,i)/domega(j,3,i),j=1,3)
 707         write (iout,*)
 708       enddo
 709       return
 710       end
 711 c------------------------------------------------------------
 712       subroutine chainbuild_cart
 713       implicit real*8 (a-h,o-z)
 714       include 'DIMENSIONS'
 715 #ifdef MPI
 716       include 'mpif.h'
 717 #endif
 718       include 'COMMON.SETUP'
 719       include 'COMMON.CHAIN'
 720       include 'COMMON.LOCAL'
 721       include 'COMMON.TIME1'
 722       include 'COMMON.IOUNITS'
 723
 724 #ifdef MPI
 725       if (nfgtasks.gt.1) then
 726 c        write (iout,*) "BCAST in chainbuild_cart"
 727 c        call flush(iout)
 728 c Broadcast the order to build the chain and compute internal coordinates
 729 c to the slaves. The slaves receive the order in ERGASTULUM.
 730         time00=MPI_Wtime()
 731 c      write (iout,*) "CHAINBUILD_CART: DC before BCAST"
 732 c      do i=0,nres
 733 c        write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
 734 c     &   (dc(j,i+nres),j=1,3)
 735 c      enddo
 736         if (fg_rank.eq.0)
 737      &    call MPI_Bcast(7,1,MPI_INTEGER,king,FG_COMM,IERROR)
 738         time_bcast7=time_bcast7+MPI_Wtime()-time00
 739         time01=MPI_Wtime()
 740         call MPI_Bcast(dc(1,0),6*(nres+1),MPI_DOUBLE_PRECISION,
 741      &    king,FG_COMM,IERR)
 742 c      write (iout,*) "CHAINBUILD_CART: DC after BCAST"
 743 c      do i=0,nres
 744 c        write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
 745 c     &   (dc(j,i+nres),j=1,3)
 746 c      enddo
 747 c        write (iout,*) "End BCAST in chainbuild_cart"
 748 c        call flush(iout)
 749         time_bcast=time_bcast+MPI_Wtime()-time00
 750         time_bcastc=time_bcastc+MPI_Wtime()-time01
 751       endif
 752 #endif
 753       do j=1,3
 754         c(j,1)=dc(j,0)
 755 c        c(j,1)=c(j,1)
 756       enddo
 757       do i=2,nres
 758         do j=1,3
 759           c(j,i)=c(j,i-1)+dc(j,i-1)
 760         enddo
 761       enddo
 762       do i=1,nres
 763         do j=1,3
 764           c(j,i+nres)=c(j,i)+dc(j,i+nres)
 765         enddo
 766       enddo
 767 C      print *,'tutu'
 768 c      write (iout,*) "CHAINBUILD_CART"
 769 c      call cartprint
 770       call int_from_cart1(.false.)
 771       return
 772       end