X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=blobdiff_plain;f=source%2Funres%2Fsrc_CSA_DiL%2Fintcartderiv.F;fp=source%2Funres%2Fsrc_CSA_DiL%2Fintcartderiv.F;h=0000000000000000000000000000000000000000;hp=5fea8759c19347c9dfca1c64982e2ef02ea1aae0;hb=2a226bfc86eabc6e4eae0c3ad1cbc3cb5417a05a;hpb=a0e685f844163003749ba91dfbf4644bcc8cfa30 diff --git a/source/unres/src_CSA_DiL/intcartderiv.F b/source/unres/src_CSA_DiL/intcartderiv.F deleted file mode 100644 index 5fea875..0000000 --- a/source/unres/src_CSA_DiL/intcartderiv.F +++ /dev/null @@ -1,466 +0,0 @@ - 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' - 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 - -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 -#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 -cd write (iout,*) "Gather dtheta" -cd call flush(iout) - write (iout,*) "dtheta before gather" - do i=1,nres - write (iout,'(i3,3(3f8.5,3x))') i,((dtheta(j,k,i),k=1,3),j=1,2) - 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