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