X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fintcartderiv.F;h=369a4f06a6476804c45f521da4cb797bde84313f;hb=2acc9917f7452c3a59a242e6033bea16f846cbfe;hp=33d4a5046a4f58100edc59d629e2df355508acf3;hpb=f5cea66c6a23f13718f902416e527a955693d213;p=unres.git diff --git a/source/unres/src_MD-M/intcartderiv.F b/source/unres/src_MD-M/intcartderiv.F index 33d4a50..369a4f0 100644 --- a/source/unres/src_MD-M/intcartderiv.F +++ b/source/unres/src_MD-M/intcartderiv.F @@ -12,6 +12,7 @@ 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), @@ -53,7 +54,42 @@ c We need dtheta(:,:,i-1) to compute dphi(:,:,i) if (itype(i-1).ne.ntyp1) dtheta(j,2,i)=-dcostheta(j,2,i)/sint 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.ntyp1)) 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 @@ -123,6 +159,228 @@ c Obtaining the gamma derivatives from cosine derivative enddo endif enddo +Calculate derivative of Tauangle +#ifdef PARINTDER + do i=itau_start,itau_end +#else + do i=3,nres +#endif + if ((itype(i-2).eq.ntyp1).or.(itype(i-2).eq.10)) cycle +c if ((itype(i-2).eq.ntyp1).or.(itype(i-2).eq.10).or. +c & (itype(i-1).eq.ntyp1).or.(itype(i).eq.ntyp1)) 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)) + 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 + 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.ntyp1).or.(itype(i-1).eq.10).or. + & (itype(i-2).eq.ntyp1).or.(itype(i-3).eq.ntyp1)) 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.ntyp1).or.(itype(i-1).eq.10).or. + &(itype(i-2).eq.ntyp1).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)