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