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