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).ge.-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+nres)=-dc_norm(j,i+nres) 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)) C do j=1,3 C dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres) cc write(iout,*) dc_norm2(j,i-2+nres),"dcnorm" C enddo 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,1,i),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",1,i,k,(dtauangle(j,1,k,i),j=1,3) c write(iout,*) "tu",2,i,k,(dtauangle(j,2,k,i),j=1,3) c write(iout,*) "tu",3,i,k,(dtauangle(j,3,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 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)) C do j=1,3 C 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) C 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