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