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