source/unres/src_MD-M/intcartderiv.F

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