--- /dev/null
+ subroutine intcartderiv
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'COMMON.SETUP'
+ include 'COMMON.CHAIN'
+ include 'COMMON.VAR'
+ include 'COMMON.GEO'
+ include 'COMMON.INTERACT'
+ include 'COMMON.DERIV'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.LOCAL'
+ include 'COMMON.SCCOR'
+ double precision dcostheta(3,2,maxres),
+ & dcosphi(3,3,maxres),dsinphi(3,3,maxres),
+ & dcosalpha(3,3,maxres),dcosomega(3,3,maxres),
+ & dsinomega(3,3,maxres),vo1(3),vo2(3),vo3(3),
+ & dummy(3),vp1(3),vp2(3),vp3(3),vpp1(3),n(3)
+
+#if defined(MPI) && defined(PARINTDER)
+ if (nfgtasks.gt.1 .and. me.eq.king)
+ & call MPI_Bcast(8,1,MPI_INTEGER,king,FG_COMM,IERROR)
+#endif
+ pi4 = 0.5d0*pipol
+ pi34 = 3*pi4
+
+c write (iout,*) "iphi1_start",iphi1_start," iphi1_end",iphi1_end
+c Derivatives of theta's
+#if defined(MPI) && defined(PARINTDER)
+c We need dtheta(:,:,i-1) to compute dphi(:,:,i)
+ do i=max0(ithet_start-1,3),ithet_end
+#else
+ do i=3,nres
+#endif
+ cost=dcos(theta(i))
+ sint=sqrt(1-cost*cost)
+ do j=1,3
+ dcostheta(j,1,i)=-(dc_norm(j,i-1)+cost*dc_norm(j,i-2))/
+ & vbld(i-1)
+ dtheta(j,1,i)=-1/sint*dcostheta(j,1,i)
+ dcostheta(j,2,i)=-(dc_norm(j,i-2)+cost*dc_norm(j,i-1))/
+ & vbld(i)
+ dtheta(j,2,i)=-1/sint*dcostheta(j,2,i)
+ enddo
+ enddo
+
+#if defined(MPI) && defined(PARINTDER)
+c We need dtheta(:,:,i-1) to compute dphi(:,:,i)
+ do i=max0(ithet_start-1,3),ithet_end
+#else
+ do i=3,nres
+#endif
+ if ((itype(i-1).ne.10).and.(itype(i-1).ne.21)) then
+ cost1=dcos(omicron(1,i))
+ sint1=sqrt(1-cost1*cost1)
+ cost2=dcos(omicron(2,i))
+ sint2=sqrt(1-cost2*cost2)
+ do j=1,3
+CC Calculate derivative over first omicron (Cai-2,Cai-1,SCi-1)
+ dcosomicron(j,1,1,i)=-(dc_norm(j,i-1+nres)+
+ & cost1*dc_norm(j,i-2))/
+ & vbld(i-1)
+ domicron(j,1,1,i)=-1/sint1*dcosomicron(j,1,1,i)
+ dcosomicron(j,1,2,i)=-(dc_norm(j,i-2)
+ & +cost1*(dc_norm(j,i-1+nres)))/
+ & vbld(i-1+nres)
+ domicron(j,1,2,i)=-1/sint1*dcosomicron(j,1,2,i)
+CC Calculate derivative over second omicron Sci-1,Cai-1 Cai
+CC Looks messy but better than if in loop
+ dcosomicron(j,2,1,i)=-(-dc_norm(j,i-1+nres)
+ & +cost2*dc_norm(j,i-1))/
+ & vbld(i)
+ domicron(j,2,1,i)=-1/sint2*dcosomicron(j,2,1,i)
+ dcosomicron(j,2,2,i)=-(dc_norm(j,i-1)
+ & +cost2*(-dc_norm(j,i-1+nres)))/
+ & vbld(i-1+nres)
+c write(iout,*) "vbld", i,itype(i),vbld(i-1+nres)
+ domicron(j,2,2,i)=-1/sint2*dcosomicron(j,2,2,i)
+ enddo
+ endif
+ enddo
+
+
+
+c Derivatives of phi:
+c If phi is 0 or 180 degrees, then the formulas
+c have to be derived by power series expansion of the
+c conventional formulas around 0 and 180.
+#ifdef PARINTDER
+ do i=iphi1_start,iphi1_end
+#else
+ do i=4,nres
+#endif
+c the conventional case
+ sint=dsin(theta(i))
+ sint1=dsin(theta(i-1))
+ sing=dsin(phi(i))
+ cost=dcos(theta(i))
+ cost1=dcos(theta(i-1))
+ cosg=dcos(phi(i))
+ scalp=scalar(dc_norm(1,i-3),dc_norm(1,i-1))
+ fac0=1.0d0/(sint1*sint)
+ fac1=cost*fac0
+ fac2=cost1*fac0
+ fac3=cosg*cost1/(sint1*sint1)
+ fac4=cosg*cost/(sint*sint)
+c Obtaining the gamma derivatives from sine derivative
+ if (phi(i).gt.-pi4.and.phi(i).le.pi4.or.
+ & phi(i).gt.pi34.and.phi(i).le.pi.or.
+ & phi(i).gt.-pi.and.phi(i).le.-pi34) then
+ call vecpr(dc_norm(1,i-1),dc_norm(1,i-2),vp1)
+ call vecpr(dc_norm(1,i-3),dc_norm(1,i-1),vp2)
+ call vecpr(dc_norm(1,i-3),dc_norm(1,i-2),vp3)
+ do j=1,3
+ ctgt=cost/sint
+ ctgt1=cost1/sint1
+ cosg_inv=1.0d0/cosg
+ dsinphi(j,1,i)=-sing*ctgt1*dtheta(j,1,i-1)
+ & -(fac0*vp1(j)+sing*dc_norm(j,i-3))*vbld_inv(i-2)
+ dphi(j,1,i)=cosg_inv*dsinphi(j,1,i)
+ dsinphi(j,2,i)=
+ & -sing*(ctgt1*dtheta(j,2,i-1)+ctgt*dtheta(j,1,i))
+ & -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
+ dphi(j,2,i)=cosg_inv*dsinphi(j,2,i)
+c Bug fixed 3/24/05 (AL)
+ dsinphi(j,3,i)=-sing*ctgt*dtheta(j,2,i)
+ & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i)
+c & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
+ dphi(j,3,i)=cosg_inv*dsinphi(j,3,i)
+ enddo
+c Obtaining the gamma derivatives from cosine derivative
+ else
+ do j=1,3
+ dcosphi(j,1,i)=fac1*dcostheta(j,1,i-1)+fac3*
+ & dcostheta(j,1,i-1)-fac0*(dc_norm(j,i-1)-scalp*
+ & dc_norm(j,i-3))/vbld(i-2)
+ dphi(j,1,i)=-1/sing*dcosphi(j,1,i)
+ dcosphi(j,2,i)=fac1*dcostheta(j,2,i-1)+fac2*
+ & dcostheta(j,1,i)+fac3*dcostheta(j,2,i-1)+fac4*
+ & dcostheta(j,1,i)
+ dphi(j,2,i)=-1/sing*dcosphi(j,2,i)
+ dcosphi(j,3,i)=fac2*dcostheta(j,2,i)+fac4*
+ & dcostheta(j,2,i)-fac0*(dc_norm(j,i-3)-scalp*
+ & dc_norm(j,i-1))/vbld(i)
+ dphi(j,3,i)=-1/sing*dcosphi(j,3,i)
+ enddo
+ endif
+ enddo
+
+ do i=1,nres-1
+ do j=1,3
+ dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres)
+cc write(iout,*) dc_norm2(j,i-2+nres),"dcnorm"
+ enddo
+ enddo
+Calculate derivative of Tauangle
+#ifdef PARINTDER
+ do i=itau_start,itau_end
+#else
+ do i=3,nres
+#endif
+ if ((itype(i-2).eq.21).or.(itype(i-2).eq.10)) cycle
+cc dtauangle(j,intertyp,dervityp,residue number)
+cc INTERTYP=1 SC...Ca...Ca..Ca
+c the conventional case
+ sint=dsin(theta(i))
+ sint1=dsin(omicron(2,i-1))
+ sing=dsin(tauangle(1,i))
+ cost=dcos(theta(i))
+ cost1=dcos(omicron(2,i-1))
+ cosg=dcos(tauangle(1,i))
+ scalp=scalar(dc_norm2(1,i-2+nres),dc_norm(1,i-1))
+ fac0=1.0d0/(sint1*sint)
+ fac1=cost*fac0
+ fac2=cost1*fac0
+ fac3=cosg*cost1/(sint1*sint1)
+ fac4=cosg*cost/(sint*sint)
+cc write(iout,*) "faki",fac0,fac1,fac2,fac3,fac4
+c Obtaining the gamma derivatives from sine derivative
+ if (tauangle(1,i).gt.-pi4.and.tauangle(1,i).le.pi4.or.
+ & tauangle(1,i).gt.pi34.and.tauangle(1,i).le.pi.or.
+ & tauangle(1,i).gt.-pi.and.tauangle(1,i).le.-pi34) then
+ call vecpr(dc_norm(1,i-1),dc_norm(1,i-2),vp1)
+ call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-1),vp2)
+ call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-2),vp3)
+ do j=1,3
+ ctgt=cost/sint
+ ctgt1=cost1/sint1
+ cosg_inv=1.0d0/cosg
+ dsintau(j,1,1,i)=-sing*ctgt1*domicron(j,2,2,i-1)
+ &-(fac0*vp1(j)+sing*(dc_norm2(j,i-2+nres)))
+ & *vbld_inv(i-2+nres)
+ dtauangle(j,1,1,i)=cosg_inv*dsintau(j,1,1,i)
+ dsintau(j,1,2,i)=
+ & -sing*(ctgt1*domicron(j,2,1,i-1)+ctgt*dtheta(j,1,i))
+ & -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
+c write(iout,*) "dsintau", dsintau(j,1,2,i)
+ dtauangle(j,1,2,i)=cosg_inv*dsintau(j,1,2,i)
+c Bug fixed 3/24/05 (AL)
+ dsintau(j,1,3,i)=-sing*ctgt*dtheta(j,2,i)
+ & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i)
+c & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
+ dtauangle(j,1,3,i)=cosg_inv*dsintau(j,1,3,i)
+ enddo
+c Obtaining the gamma derivatives from cosine derivative
+ else
+ do j=1,3
+ dcostau(j,1,1,i)=fac1*dcosomicron(j,2,2,i-1)+fac3*
+ & dcosomicron(j,2,2,i-1)-fac0*(dc_norm(j,i-1)-scalp*
+ & (dc_norm2(j,i-2+nres)))/vbld(i-2+nres)
+ dtauangle(j,1,1,i)=-1/sing*dcostau(j,1,1,i)
+ dcostau(j,1,2,i)=fac1*dcosomicron(j,2,1,i-1)+fac2*
+ & dcostheta(j,1,i)+fac3*dcosomicron(j,2,1,i-1)+fac4*
+ & dcostheta(j,1,i)
+ dtauangle(j,1,2,i)=-1/sing*dcostau(j,1,2,i)
+ dcostau(j,1,3,i)=fac2*dcostheta(j,2,i)+fac4*
+ & dcostheta(j,2,i)-fac0*(-dc_norm(j,i-2+nres)-scalp*
+ & dc_norm(j,i-1))/vbld(i)
+ dtauangle(j,1,3,i)=-1/sing*dcostau(j,1,3,i)
+c write (iout,*) "else",i
+ enddo
+ endif
+c do k=1,3
+c write(iout,*) "tu",i,k,(dtauangle(j,1,k,i),j=1,3)
+c enddo
+ enddo
+CC Second case Ca...Ca...Ca...SC
+#ifdef PARINTDER
+ do i=itau_start,itau_end
+#else
+ do i=4,nres
+#endif
+ if ((itype(i-1).eq.21).or.(itype(i-1).eq.10)) cycle
+c the conventional case
+ sint=dsin(omicron(1,i))
+ sint1=dsin(theta(i-1))
+ sing=dsin(tauangle(2,i))
+ cost=dcos(omicron(1,i))
+ cost1=dcos(theta(i-1))
+ cosg=dcos(tauangle(2,i))
+c do j=1,3
+c dc_norm2(j,i-1+nres)=-dc_norm(j,i-1+nres)
+c enddo
+ scalp=scalar(dc_norm(1,i-3),dc_norm(1,i-1+nres))
+ fac0=1.0d0/(sint1*sint)
+ fac1=cost*fac0
+ fac2=cost1*fac0
+ fac3=cosg*cost1/(sint1*sint1)
+ fac4=cosg*cost/(sint*sint)
+c Obtaining the gamma derivatives from sine derivative
+c write (iout,*) "i",i," tauangle2",tauangle(2,i)
+ if (tauangle(2,i).gt.-pi4.and.tauangle(2,i).le.pi4.or.
+ & tauangle(2,i).gt.pi34.and.tauangle(2,i).le.pi.or.
+ & tauangle(2,i).gt.-pi.and.tauangle(2,i).le.-pi34) then
+ call vecpr(dc_norm2(1,i-1+nres),dc_norm(1,i-2),vp1)
+ call vecpr(dc_norm(1,i-3),dc_norm(1,i-1+nres),vp2)
+ call vecpr(dc_norm(1,i-3),dc_norm(1,i-2),vp3)
+ do j=1,3
+ ctgt=cost/sint
+ ctgt1=cost1/sint1
+ cosg_inv=1.0d0/cosg
+ dsintau(j,2,1,i)=-sing*ctgt1*dtheta(j,1,i-1)
+ & +(fac0*vp1(j)-sing*dc_norm(j,i-3))*vbld_inv(i-2)
+c write(iout,*) i,j,dsintau(j,2,1,i),sing*ctgt1*dtheta(j,1,i-1),
+c &fac0*vp1(j),sing*dc_norm(j,i-3),vbld_inv(i-2),"dsintau(2,1)"
+ dtauangle(j,2,1,i)=cosg_inv*dsintau(j,2,1,i)
+ dsintau(j,2,2,i)=
+ & -sing*(ctgt1*dtheta(j,2,i-1)+ctgt*domicron(j,1,1,i))
+ & -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
+c write(iout,*) "sprawdzenie",i,j,sing*ctgt1*dtheta(j,2,i-1),
+c & sing*ctgt*domicron(j,1,2,i),
+c & (fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
+ dtauangle(j,2,2,i)=cosg_inv*dsintau(j,2,2,i)
+c Bug fixed 3/24/05 (AL)
+ dsintau(j,2,3,i)=-sing*ctgt*domicron(j,1,2,i)
+ & +(fac0*vp3(j)-sing*dc_norm(j,i-1+nres))*vbld_inv(i-1+nres)
+c & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
+ dtauangle(j,2,3,i)=cosg_inv*dsintau(j,2,3,i)
+ enddo
+c Obtaining the gamma derivatives from cosine derivative
+ else
+ do j=1,3
+ dcostau(j,2,1,i)=fac1*dcostheta(j,1,i-1)+fac3*
+ & dcostheta(j,1,i-1)-fac0*(dc_norm(j,i-1+nres)-scalp*
+ & dc_norm(j,i-3))/vbld(i-2)
+ dtauangle(j,2,1,i)=-1/sing*dcostau(j,2,1,i)
+ dcostau(j,2,2,i)=fac1*dcostheta(j,2,i-1)+fac2*
+ & dcosomicron(j,1,1,i)+fac3*dcostheta(j,2,i-1)+fac4*
+ & dcosomicron(j,1,1,i)
+ dtauangle(j,2,2,i)=-1/sing*dcostau(j,2,2,i)
+ dcostau(j,2,3,i)=fac2*dcosomicron(j,1,2,i)+fac4*
+ & dcosomicron(j,1,2,i)-fac0*(dc_norm(j,i-3)-scalp*
+ & dc_norm(j,i-1+nres))/vbld(i-1+nres)
+ dtauangle(j,2,3,i)=-1/sing*dcostau(j,2,3,i)
+c write(iout,*) i,j,"else", dtauangle(j,2,3,i)
+ enddo
+ endif
+ enddo
+
+
+CCC third case SC...Ca...Ca...SC
+#ifdef PARINTDER
+
+ do i=itau_start,itau_end
+#else
+ do i=3,nres
+#endif
+c the conventional case
+ if ((itype(i-1).eq.21).or.(itype(i-1).eq.10).or.
+ &(itype(i-2).eq.21).or.(itype(i-2).eq.10)) cycle
+ sint=dsin(omicron(1,i))
+ sint1=dsin(omicron(2,i-1))
+ sing=dsin(tauangle(3,i))
+ cost=dcos(omicron(1,i))
+ cost1=dcos(omicron(2,i-1))
+ cosg=dcos(tauangle(3,i))
+ do j=1,3
+ dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres)
+c dc_norm2(j,i-1+nres)=-dc_norm(j,i-1+nres)
+ enddo
+ scalp=scalar(dc_norm2(1,i-2+nres),dc_norm(1,i-1+nres))
+ fac0=1.0d0/(sint1*sint)
+ fac1=cost*fac0
+ fac2=cost1*fac0
+ fac3=cosg*cost1/(sint1*sint1)
+ fac4=cosg*cost/(sint*sint)
+c Obtaining the gamma derivatives from sine derivative
+ if (tauangle(3,i).gt.-pi4.and.tauangle(3,i).le.pi4.or.
+ & tauangle(3,i).gt.pi34.and.tauangle(3,i).le.pi.or.
+ & tauangle(3,i).gt.-pi.and.tauangle(3,i).le.-pi34) then
+ call vecpr(dc_norm(1,i-1+nres),dc_norm(1,i-2),vp1)
+ call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-1+nres),vp2)
+ call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-2),vp3)
+ do j=1,3
+ ctgt=cost/sint
+ ctgt1=cost1/sint1
+ cosg_inv=1.0d0/cosg
+ dsintau(j,3,1,i)=-sing*ctgt1*domicron(j,2,2,i-1)
+ & -(fac0*vp1(j)-sing*dc_norm(j,i-2+nres))
+ & *vbld_inv(i-2+nres)
+ dtauangle(j,3,1,i)=cosg_inv*dsintau(j,3,1,i)
+ dsintau(j,3,2,i)=
+ & -sing*(ctgt1*domicron(j,2,1,i-1)+ctgt*domicron(j,1,1,i))
+ & -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
+ dtauangle(j,3,2,i)=cosg_inv*dsintau(j,3,2,i)
+c Bug fixed 3/24/05 (AL)
+ dsintau(j,3,3,i)=-sing*ctgt*domicron(j,1,2,i)
+ & +(fac0*vp3(j)-sing*dc_norm(j,i-1+nres))
+ & *vbld_inv(i-1+nres)
+c & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
+ dtauangle(j,3,3,i)=cosg_inv*dsintau(j,3,3,i)
+ enddo
+c Obtaining the gamma derivatives from cosine derivative
+ else
+ do j=1,3
+ dcostau(j,3,1,i)=fac1*dcosomicron(j,2,2,i-1)+fac3*
+ & dcosomicron(j,2,2,i-1)-fac0*(dc_norm(j,i-1+nres)-scalp*
+ & dc_norm2(j,i-2+nres))/vbld(i-2+nres)
+ dtauangle(j,3,1,i)=-1/sing*dcostau(j,3,1,i)
+ dcostau(j,3,2,i)=fac1*dcosomicron(j,2,1,i-1)+fac2*
+ & dcosomicron(j,1,1,i)+fac3*dcosomicron(j,2,1,i-1)+fac4*
+ & dcosomicron(j,1,1,i)
+ dtauangle(j,3,2,i)=-1/sing*dcostau(j,3,2,i)
+ dcostau(j,3,3,i)=fac2*dcosomicron(j,1,2,i)+fac4*
+ & dcosomicron(j,1,2,i)-fac0*(dc_norm2(j,i-2+nres)-scalp*
+ & dc_norm(j,i-1+nres))/vbld(i-1+nres)
+ dtauangle(j,3,3,i)=-1/sing*dcostau(j,3,3,i)
+c write(iout,*) "else",i
+ enddo
+ endif
+ enddo
+#ifdef CRYST_SC
+c Derivatives of side-chain angles alpha and omega
+#if defined(MPI) && defined(PARINTDER)
+ do i=ibond_start,ibond_end
+#else
+ do i=2,nres-1
+#endif
+ if(itype(i).ne.10) then
+ fac5=1.0d0/dsqrt(2*(1+dcos(theta(i+1))))
+ fac6=fac5/vbld(i)
+ fac7=fac5*fac5
+ fac8=fac5/vbld(i+1)
+ fac9=fac5/vbld(i+nres)
+ scala1=scalar(dc_norm(1,i-1),dc_norm(1,i+nres))
+ scala2=scalar(dc_norm(1,i),dc_norm(1,i+nres))
+ cosa=dsqrt(0.5d0/(1.0d0+dcos(theta(i+1))))*(
+ & scalar(dC_norm(1,i),dC_norm(1,i+nres))
+ & -scalar(dC_norm(1,i-1),dC_norm(1,i+nres)))
+ sina=sqrt(1-cosa*cosa)
+ sino=dsin(omeg(i))
+ do j=1,3
+ dcosalpha(j,1,i)=fac6*(scala1*dc_norm(j,i-1)-
+ & dc_norm(j,i+nres))-cosa*fac7*dcostheta(j,1,i+1)
+ dalpha(j,1,i)=-1/sina*dcosalpha(j,1,i)
+ dcosalpha(j,2,i)=fac8*(dc_norm(j,i+nres)-
+ & scala2*dc_norm(j,i))-cosa*fac7*dcostheta(j,2,i+1)
+ dalpha(j,2,i)=-1/sina*dcosalpha(j,2,i)
+ dcosalpha(j,3,i)=(fac9*(dc_norm(j,i)-
+ & dc_norm(j,i-1))-(cosa*dc_norm(j,i+nres))/
+ & vbld(i+nres))
+ dalpha(j,3,i)=-1/sina*dcosalpha(j,3,i)
+ enddo
+c obtaining the derivatives of omega from sines
+ if(omeg(i).gt.-pi4.and.omeg(i).le.pi4.or.
+ & omeg(i).gt.pi34.and.omeg(i).le.pi.or.
+ & omeg(i).gt.-pi.and.omeg(i).le.-pi34) then
+ fac15=dcos(theta(i+1))/(dsin(theta(i+1))*
+ & dsin(theta(i+1)))
+ fac16=dcos(alph(i))/(dsin(alph(i))*dsin(alph(i)))
+ fac17=1.0d0/(dsin(theta(i+1))*dsin(alph(i)))
+ call vecpr(dc_norm(1,i+nres),dc_norm(1,i),vo1)
+ call vecpr(dc_norm(1,i+nres),dc_norm(1,i-1),vo2)
+ call vecpr(dc_norm(1,i),dc_norm(1,i-1),vo3)
+ coso_inv=1.0d0/dcos(omeg(i))
+ do j=1,3
+ dsinomega(j,1,i)=sino*(fac15*dcostheta(j,1,i+1)
+ & +fac16*dcosalpha(j,1,i))-fac17/vbld(i)*vo1(j)-(
+ & sino*dc_norm(j,i-1))/vbld(i)
+ domega(j,1,i)=coso_inv*dsinomega(j,1,i)
+ dsinomega(j,2,i)=sino*(fac15*dcostheta(j,2,i+1)
+ & +fac16*dcosalpha(j,2,i))+fac17/vbld(i+1)*vo2(j)
+ & -sino*dc_norm(j,i)/vbld(i+1)
+ domega(j,2,i)=coso_inv*dsinomega(j,2,i)
+ dsinomega(j,3,i)=sino*fac16*dcosalpha(j,3,i)-
+ & fac17/vbld(i+nres)*vo3(j)-sino*dc_norm(j,i+nres)/
+ & vbld(i+nres)
+ domega(j,3,i)=coso_inv*dsinomega(j,3,i)
+ enddo
+ else
+c obtaining the derivatives of omega from cosines
+ fac10=sqrt(0.5d0*(1-dcos(theta(i+1))))
+ fac11=sqrt(0.5d0*(1+dcos(theta(i+1))))
+ fac12=fac10*sina
+ fac13=fac12*fac12
+ fac14=sina*sina
+ do j=1,3
+ dcosomega(j,1,i)=(-(0.25d0*cosa/fac11*
+ & dcostheta(j,1,i+1)+fac11*dcosalpha(j,1,i))*fac12+
+ & (0.25d0/fac10*sina*dcostheta(j,1,i+1)+cosa/sina*
+ & fac10*dcosalpha(j,1,i))*(scala2-fac11*cosa))/fac13
+ domega(j,1,i)=-1/sino*dcosomega(j,1,i)
+ dcosomega(j,2,i)=(((dc_norm(j,i+nres)-scala2*
+ & dc_norm(j,i))/vbld(i+1)-0.25d0*cosa/fac11*
+ & dcostheta(j,2,i+1)-fac11*dcosalpha(j,2,i))*fac12+
+ & (scala2-fac11*cosa)*(0.25d0*sina/fac10*
+ & dcostheta(j,2,i+1)+fac10*cosa/sina*dcosalpha(j,2,i)
+ & ))/fac13
+ domega(j,2,i)=-1/sino*dcosomega(j,2,i)
+ dcosomega(j,3,i)=1/fac10*((1/vbld(i+nres)*(dc_norm(j,i)-
+ & scala2*dc_norm(j,i+nres))-fac11*dcosalpha(j,3,i))*sina+
+ & (scala2-fac11*cosa)*(cosa/sina*dcosalpha(j,3,i)))/fac14
+ domega(j,3,i)=-1/sino*dcosomega(j,3,i)
+ enddo
+ endif
+ endif
+ enddo
+#endif
+#if defined(MPI) && defined(PARINTDER)
+ if (nfgtasks.gt.1) then
+#ifdef DEBUG
+ write (iout,*) "Gather dtheta"
+cd call flush(iout)
+c write (iout,*) "dtheta before gather"
+c do i=1,nres
+c write (iout,'(i3,3(3f8.5,3x))') i,((dtheta(j,k,i),k=1,3),j=1,2)
+c enddo
+#endif
+ call MPI_Gatherv(dtheta(1,1,ithet_start),ithet_count(fg_rank),
+ & MPI_THET,dtheta(1,1,1),ithet_count(0),ithet_displ(0),MPI_THET,
+ & king,FG_COMM,IERROR)
+#ifdef DEBUG
+cd write (iout,*) "Gather dphi"
+cd call flush(iout)
+ write (iout,*) "dphi before gather"
+ do i=1,nres
+ write (iout,'(i3,3(3f8.5,3x))') i,((dphi(j,k,i),k=1,3),j=1,3)
+ enddo
+#endif
+ call MPI_Gatherv(dphi(1,1,iphi1_start),iphi1_count(fg_rank),
+ & MPI_GAM,dphi(1,1,1),iphi1_count(0),iphi1_displ(0),MPI_GAM,
+ & king,FG_COMM,IERROR)
+cd write (iout,*) "Gather dalpha"
+cd call flush(iout)
+#ifdef CRYST_SC
+ call MPI_Gatherv(dalpha(1,1,ibond_start),ibond_count(fg_rank),
+ & MPI_GAM,dalpha(1,1,1),ibond_count(0),ibond_displ(0),MPI_GAM,
+ & king,FG_COMM,IERROR)
+cd write (iout,*) "Gather domega"
+cd call flush(iout)
+ call MPI_Gatherv(domega(1,1,ibond_start),ibond_count(fg_rank),
+ & MPI_GAM,domega(1,1,1),ibond_count(0),ibond_displ(0),MPI_GAM,
+ & king,FG_COMM,IERROR)
+#endif
+ endif
+#endif
+#ifdef DEBUG
+ write (iout,*) "dtheta after gather"
+ do i=1,nres
+ write (iout,'(i3,3(3f8.5,3x))') i,((dtheta(j,k,i),j=1,3),j=1,2)
+ enddo
+ write (iout,*) "dphi after gather"
+ do i=1,nres
+ write (iout,'(i3,3(3f8.5,3x))') i,((dphi(j,k,i),j=1,3),k=1,3)
+ enddo
+#endif
+ return
+ end
+
+ subroutine checkintcartgrad
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'COMMON.CHAIN'
+ include 'COMMON.VAR'
+ include 'COMMON.GEO'
+ include 'COMMON.INTERACT'
+ include 'COMMON.DERIV'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.SETUP'
+ double precision dthetanum(3,2,maxres),dphinum(3,3,maxres)
+ & ,dalphanum(3,3,maxres), domeganum(3,3,maxres)
+ double precision theta_s(maxres),phi_s(maxres),alph_s(maxres),
+ & omeg_s(maxres),dc_norm_s(3)
+ double precision aincr /1.0d-5/
+
+ do i=1,nres
+ phi_s(i)=phi(i)
+ theta_s(i)=theta(i)
+ alph_s(i)=alph(i)
+ omeg_s(i)=omeg(i)
+ enddo
+c Check theta gradient
+ write (iout,*)
+ & "Analytical (upper) and numerical (lower) gradient of theta"
+ write (iout,*)
+ do i=3,nres
+ do j=1,3
+ dcji=dc(j,i-2)
+ dc(j,i-2)=dcji+aincr
+ call chainbuild_cart
+ call int_from_cart1(.false.)
+ dthetanum(j,1,i)=(theta(i)-theta_s(i))/aincr
+ dc(j,i-2)=dcji
+ dcji=dc(j,i-1)
+ dc(j,i-1)=dc(j,i-1)+aincr
+ call chainbuild_cart
+ dthetanum(j,2,i)=(theta(i)-theta_s(i))/aincr
+ dc(j,i-1)=dcji
+ enddo
+ write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dtheta(j,1,i),j=1,3),
+ & (dtheta(j,2,i),j=1,3)
+ write (iout,'(5x,3f10.5,5x,3f10.5)') (dthetanum(j,1,i),j=1,3),
+ & (dthetanum(j,2,i),j=1,3)
+ write (iout,'(5x,3f10.5,5x,3f10.5)')
+ & (dthetanum(j,1,i)/dtheta(j,1,i),j=1,3),
+ & (dthetanum(j,2,i)/dtheta(j,2,i),j=1,3)
+ write (iout,*)
+ enddo
+c Check gamma gradient
+ write (iout,*)
+ & "Analytical (upper) and numerical (lower) gradient of gamma"
+ do i=4,nres
+ do j=1,3
+ dcji=dc(j,i-3)
+ dc(j,i-3)=dcji+aincr
+ call chainbuild_cart
+ dphinum(j,1,i)=(phi(i)-phi_s(i))/aincr
+ dc(j,i-3)=dcji
+ dcji=dc(j,i-2)
+ dc(j,i-2)=dcji+aincr
+ call chainbuild_cart
+ dphinum(j,2,i)=(phi(i)-phi_s(i))/aincr
+ dc(j,i-2)=dcji
+ dcji=dc(j,i-1)
+ dc(j,i-1)=dc(j,i-1)+aincr
+ call chainbuild_cart
+ dphinum(j,3,i)=(phi(i)-phi_s(i))/aincr
+ dc(j,i-1)=dcji
+ enddo
+ write (iout,'(i5,3(3f10.5,5x))') i,(dphi(j,1,i),j=1,3),
+ & (dphi(j,2,i),j=1,3),(dphi(j,3,i),j=1,3)
+ write (iout,'(5x,3(3f10.5,5x))') (dphinum(j,1,i),j=1,3),
+ & (dphinum(j,2,i),j=1,3),(dphinum(j,3,i),j=1,3)
+ write (iout,'(5x,3(3f10.5,5x))')
+ & (dphinum(j,1,i)/dphi(j,1,i),j=1,3),
+ & (dphinum(j,2,i)/dphi(j,2,i),j=1,3),
+ & (dphinum(j,3,i)/dphi(j,3,i),j=1,3)
+ write (iout,*)
+ enddo
+c Check alpha gradient
+ write (iout,*)
+ & "Analytical (upper) and numerical (lower) gradient of alpha"
+ do i=2,nres-1
+ if(itype(i).ne.10) then
+ do j=1,3
+ dcji=dc(j,i-1)
+ dc(j,i-1)=dcji+aincr
+ call chainbuild_cart
+ dalphanum(j,1,i)=(alph(i)-alph_s(i))
+ & /aincr
+ dc(j,i-1)=dcji
+ dcji=dc(j,i)
+ dc(j,i)=dcji+aincr
+ call chainbuild_cart
+ dalphanum(j,2,i)=(alph(i)-alph_s(i))
+ & /aincr
+ dc(j,i)=dcji
+ dcji=dc(j,i+nres)
+ dc(j,i+nres)=dc(j,i+nres)+aincr
+ call chainbuild_cart
+ dalphanum(j,3,i)=(alph(i)-alph_s(i))
+ & /aincr
+ dc(j,i+nres)=dcji
+ enddo
+ endif
+ write (iout,'(i5,3(3f10.5,5x))') i,(dalpha(j,1,i),j=1,3),
+ & (dalpha(j,2,i),j=1,3),(dalpha(j,3,i),j=1,3)
+ write (iout,'(5x,3(3f10.5,5x))') (dalphanum(j,1,i),j=1,3),
+ & (dalphanum(j,2,i),j=1,3),(dalphanum(j,3,i),j=1,3)
+ write (iout,'(5x,3(3f10.5,5x))')
+ & (dalphanum(j,1,i)/dalpha(j,1,i),j=1,3),
+ & (dalphanum(j,2,i)/dalpha(j,2,i),j=1,3),
+ & (dalphanum(j,3,i)/dalpha(j,3,i),j=1,3)
+ write (iout,*)
+ enddo
+c Check omega gradient
+ write (iout,*)
+ & "Analytical (upper) and numerical (lower) gradient of omega"
+ do i=2,nres-1
+ if(itype(i).ne.10) then
+ do j=1,3
+ dcji=dc(j,i-1)
+ dc(j,i-1)=dcji+aincr
+ call chainbuild_cart
+ domeganum(j,1,i)=(omeg(i)-omeg_s(i))
+ & /aincr
+ dc(j,i-1)=dcji
+ dcji=dc(j,i)
+ dc(j,i)=dcji+aincr
+ call chainbuild_cart
+ domeganum(j,2,i)=(omeg(i)-omeg_s(i))
+ & /aincr
+ dc(j,i)=dcji
+ dcji=dc(j,i+nres)
+ dc(j,i+nres)=dc(j,i+nres)+aincr
+ call chainbuild_cart
+ domeganum(j,3,i)=(omeg(i)-omeg_s(i))
+ & /aincr
+ dc(j,i+nres)=dcji
+ enddo
+ endif
+ write (iout,'(i5,3(3f10.5,5x))') i,(domega(j,1,i),j=1,3),
+ & (domega(j,2,i),j=1,3),(domega(j,3,i),j=1,3)
+ write (iout,'(5x,3(3f10.5,5x))') (domeganum(j,1,i),j=1,3),
+ & (domeganum(j,2,i),j=1,3),(domeganum(j,3,i),j=1,3)
+ write (iout,'(5x,3(3f10.5,5x))')
+ & (domeganum(j,1,i)/domega(j,1,i),j=1,3),
+ & (domeganum(j,2,i)/domega(j,2,i),j=1,3),
+ & (domeganum(j,3,i)/domega(j,3,i),j=1,3)
+ write (iout,*)
+ enddo
+ return
+ end
+
+ subroutine chainbuild_cart
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'COMMON.SETUP'
+ include 'COMMON.CHAIN'
+ include 'COMMON.LOCAL'
+ include 'COMMON.TIME1'
+ include 'COMMON.IOUNITS'
+
+#ifdef MPI
+ if (nfgtasks.gt.1) then
+c write (iout,*) "BCAST in chainbuild_cart"
+c call flush(iout)
+c Broadcast the order to build the chain and compute internal coordinates
+c to the slaves. The slaves receive the order in ERGASTULUM.
+ time00=MPI_Wtime()
+c write (iout,*) "CHAINBUILD_CART: DC before BCAST"
+c do i=0,nres
+c write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
+c & (dc(j,i+nres),j=1,3)
+c enddo
+ if (fg_rank.eq.0)
+ & call MPI_Bcast(7,1,MPI_INTEGER,king,FG_COMM,IERROR)
+ time_bcast7=time_bcast7+MPI_Wtime()-time00
+ time01=MPI_Wtime()
+ call MPI_Bcast(dc(1,0),6*(nres+1),MPI_DOUBLE_PRECISION,
+ & king,FG_COMM,IERR)
+c write (iout,*) "CHAINBUILD_CART: DC after BCAST"
+c do i=0,nres
+c write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
+c & (dc(j,i+nres),j=1,3)
+c enddo
+c write (iout,*) "End BCAST in chainbuild_cart"
+c call flush(iout)
+ time_bcast=time_bcast+MPI_Wtime()-time00
+ time_bcastc=time_bcastc+MPI_Wtime()-time01
+ endif
+#endif
+ do j=1,3
+ c(j,1)=dc(j,0)
+ enddo
+ do i=2,nres
+ do j=1,3
+ c(j,i)=c(j,i-1)+dc(j,i-1)
+ enddo
+ enddo
+ do i=1,nres
+ do j=1,3
+ c(j,i+nres)=c(j,i)+dc(j,i+nres)
+ enddo
+ enddo
+c write (iout,*) "CHAINBUILD_CART"
+c call cartprint
+ call int_from_cart1(.false.)
+ return
+ end
projects / unres.git / commitdiff