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