source/maxlik/src-Fmatch/intcartderiv.F

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