X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fgeometry.f90;h=52c3b833785307229bd2e74bc7425d4f49c7abb5;hb=8ca97b16fe25b7053f258263899ba030572cc58f;hp=8542f52c9d48888f31e07e6483cbeb8378c42b7f;hpb=7c0faf2ccb6f94ed9a77aa527d3885ba054d3fb2;p=unres4.git diff --git a/source/unres/geometry.f90 b/source/unres/geometry.f90 index 8542f52..52c3b83 100644 --- a/source/unres/geometry.f90 +++ b/source/unres/geometry.f90 @@ -450,6 +450,7 @@ alph(i)=alpha(nres+i,i,nres2+2) theta(i+1)=alpha(i-1,i,i+1) vbld(i)=dist(i-1,i) +! print *,i,vbld(i),"vbld(i)" vbld_inv(i)=1.0d0/vbld(i) vbld(nres+i)=dist(nres+i,i) if (itype(i,1).ne.10) then @@ -1617,6 +1618,68 @@ end subroutine sc_angular !----------------------------------------------------------------------------- ! initialize_p.F + subroutine sc_angular_nucl +! Calculate eps1,eps2,eps3,sigma, and parts of their derivatives in om1,om2, +! om12. Called by ebp, egb, and egbv. +! use calc_data +! implicit none +! include 'COMMON.CALC' +! include 'COMMON.IOUNITS' + use comm_locel + use calc_data_nucl + erij(1)=xj*rij + erij(2)=yj*rij + erij(3)=zj*rij + om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3) + om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3) + om12=dxi*dxj+dyi*dyj+dzi*dzj + chiom12=chi12*om12 +! Calculate eps1(om12) and its derivative in om12 + faceps1=1.0D0-om12*chiom12 + faceps1_inv=1.0D0/faceps1 + eps1=dsqrt(faceps1_inv) +! Following variable is eps1*deps1/dom12 + eps1_om12=faceps1_inv*chiom12 +! diagnostics only +! faceps1_inv=om12 +! eps1=om12 +! eps1_om12=1.0d0 +! write (iout,*) "om12",om12," eps1",eps1 +! Calculate sigma(om1,om2,om12) and the derivatives of sigma**2 in om1,om2, +! and om12. + om1om2=om1*om2 + chiom1=chi1*om1 + chiom2=chi2*om2 + facsig=om1*chiom1+om2*chiom2-2.0D0*om1om2*chiom12 + sigsq=1.0D0-facsig*faceps1_inv + sigsq_om1=(chiom1-chiom12*om2)*faceps1_inv + sigsq_om2=(chiom2-chiom12*om1)*faceps1_inv + sigsq_om12=-chi12*(om1om2*faceps1-om12*facsig)*faceps1_inv**2 + chipom1=chip1*om1 + chipom2=chip2*om2 + chipom12=chip12*om12 + facp=1.0D0-om12*chipom12 + facp_inv=1.0D0/facp + facp1=om1*chipom1+om2*chipom2-2.0D0*om1om2*chipom12 +! write (iout,*) "chipom1",chipom1," chipom2",chipom2, +! & " chipom12",chipom12," facp",facp," facp_inv",facp_inv +! Following variable is the square root of eps2 + eps2rt=1.0D0-facp1*facp_inv +! Following three variables are the derivatives of the square root of eps +! in om1, om2, and om12. + eps2rt_om1=-4.0D0*(chipom1-chipom12*om2)*facp_inv + eps2rt_om2=-4.0D0*(chipom2-chipom12*om1)*facp_inv + eps2rt_om12=4.0D0*chip12*(om1om2*facp-om12*facp1)*facp_inv**2 +! Evaluate the "asymmetric" factor in the VDW constant, eps3 + eps3rt=1.0D0-alf1*om1+alf2*om2-alf12*om12 +! write (iout,*) "eps2rt",eps2rt," eps3rt",eps3rt +! write (iout,*) "eps2rt_om1",eps2rt_om1," eps2rt_om2",eps2rt_om2, +! & " eps2rt_om12",eps2rt_om12 +! Calculate whole angle-dependent part of epsilon and contributions +! to its derivatives + return + end subroutine sc_angular_nucl + !----------------------------------------------------------------------------- subroutine int_bounds(total_ints,lower_bound,upper_bound) ! implicit real*8 (a-h,o-z) @@ -2851,10 +2914,11 @@ endif endif do i=1,nres-1 +! if (molnum(i).ne.1) cycle !in wham do i=1,nres iti=itype(i,1) - if ((dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0).and.& - (iti.ne.ntyp1 .and. itype(i+1,1).ne.ntyp1)) then + if (((dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0).and.& + (iti.ne.ntyp1 .and. itype(i+1,1).ne.ntyp1)).and.molnum(i).eq.1) then write (iout,'(a,i4)') 'Bad Cartesians for residue',i !test stop endif @@ -2903,8 +2967,9 @@ di=dist(i,nres+i) !#ifndef WHAM_RUN ! 10/03/12 Adam: Correction for zero SC-SC bond length + if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1 .and. di.eq.0.0d0) & - di=dsc(itype(i,1)) + di=dsc(itype(i,molnum(i))) vbld(i+nres)=di if (itype(i,1).ne.10) then vbld_inv(i+nres)=1.0d0/di @@ -2916,12 +2981,21 @@ alph(i)=alpha(nres+i,i,nres2+2) omeg(i)=beta(nres+i,i,nres2+2,i+1) endif + if (iti.ne.0) then if(me.eq.king.or..not.out1file)then if (lprn) & write (iout,'(a3,i4,7f10.3)') restyp(iti,1),i,vbld(i),& rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),& rad2deg*alph(i),rad2deg*omeg(i) endif + else + if(me.eq.king.or..not.out1file)then + if (lprn) & + write (iout,'(a3,i4,7f10.3)') restyp(iti,1),i,vbld(i),& + rad2deg*theta(i),rad2deg*phi(i),dsc(iti+1),vbld(nres+i),& + rad2deg*alph(i),rad2deg*omeg(i) + endif + endif enddo else if (lprn) then do i=2,nres @@ -3043,7 +3117,8 @@ do j=1,3 sccmj=0.0D0 do i=1,nscat - sccmj=sccmj+sccor(j,i) + sccmj=sccmj+sccor(j,i) +!C print *,"insccent", ires,sccor(j,i) enddo dc(j,ires)=sccmj/nscat enddo @@ -3166,11 +3241,14 @@ ! calculating dE/ddc1 !el local variables integer :: j,i +! print *,"gloc",gloc(:,:) +! print *, "gcart",gcart(:,:) if (nres.lt.3) go to 18 do j=1,3 gcart(j,1)=gcart(j,1)+gloc(1,icg)*dphi(j,1,4) & +gloc(nres-2,icg)*dtheta(j,1,3) - if(itype(2,1).ne.10) then + if ((itype(2,1).ne.10).and.& + (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)))) then gcart(j,1)=gcart(j,1)+gloc(ialph(2,1),icg)*dalpha(j,1,2)+ & gloc(ialph(2,1)+nside,icg)*domega(j,1,2) endif @@ -3197,11 +3275,15 @@ gcart(j,3)=gcart(j,3)+gloc(1,icg)*dphi(j,3,4)+gloc(2,icg)* & dphi(j,2,5)+gloc(nres-1,icg)*dtheta(j,2,4)+gloc(nres,icg)* & dtheta(j,1,5) - if(itype(3,1).ne.10) then +! if(itype(3,1).ne.10) then + if ((itype(3,1).ne.10).and.& + (itype(3,molnum(3)).ne.ntyp1_molec(molnum(3)))) then gcart(j,3)=gcart(j,3)+gloc(ialph(3,1),icg)* & dalpha(j,2,3)+gloc(ialph(3,1)+nside,icg)*domega(j,2,3) endif - if(itype(4,1).ne.10) then +! if(itype(4,1).ne.10) then + if ((itype(4,1).ne.10).and.& + (itype(4,molnum(4)).ne.ntyp1_molec(molnum(4)))) then gcart(j,3)=gcart(j,3)+gloc(ialph(4,1),icg)* & dalpha(j,1,4)+gloc(ialph(4,1)+nside,icg)*domega(j,1,4) endif @@ -3257,7 +3339,8 @@ enddo ! The side-chain vector derivatives do i=2,nres-1 - if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then + if(itype(i,1).ne.10 .and. & + itype(i,molnum(i)).ne.ntyp1_molec(molnum(i))) then do j=1,3 gxcart(j,i)=gxcart(j,i)+gloc(ialph(i,1),icg)*dalpha(j,3,i) & +gloc(ialph(i,1)+nside,icg)*domega(j,3,i) @@ -3276,7 +3359,8 @@ ! enddo if (nres.lt.2) return if ((nres.lt.3).and.(itype(1,1).eq.10)) return - if ((itype(1,1).ne.10).and.(itype(1,1).ne.ntyp1)) then + if ((itype(1,1).ne.10).and. & + (itype(1,molnum(1)).ne.ntyp1_molec(molnum(1)))) then do j=1,3 !c Derviative was calculated for oposite vector of side chain therefore ! there is "-" sign before gloc_sc @@ -3284,7 +3368,8 @@ dtauangle(j,1,1,3) gcart(j,1)=gcart(j,1)+gloc_sc(1,0,icg)* & dtauangle(j,1,2,3) - if ((itype(2,1).ne.10).and.(itype(2,1).ne.ntyp1)) then + if ((itype(2,1).ne.10).and. & + (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)))) then gxcart(j,1)= gxcart(j,1) & -gloc_sc(3,0,icg)*dtauangle(j,3,1,3) gcart(j,1)=gcart(j,1)+gloc_sc(3,0,icg)* & @@ -3292,7 +3377,8 @@ endif enddo endif - if ((nres.ge.3).and.(itype(3,1).ne.10).and.(itype(3,1).ne.ntyp1)) & + if ((nres.ge.3).and.(itype(3,molnum(3)).ne.10).and.& + (itype(3,molnum(3)).ne.ntyp1_molec(molnum(3)))) & then do j=1,3 gcart(j,1)=gcart(j,1)+gloc_sc(2,1,icg)*dtauangle(j,2,1,4) @@ -3304,16 +3390,21 @@ ! Calculating the remainder of dE/ddc2 do j=1,3 - if((itype(2,1).ne.10).and.(itype(2,1).ne.ntyp1)) then - if (itype(1,1).ne.10) gxcart(j,2)=gxcart(j,2)+ & + if((itype(2,1).ne.10).and. & + (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)))) then + if ((itype(1,1).ne.10).and.& + ((itype(1,molnum(1)).ne.ntyp1_molec(molnum(1)))))& + gxcart(j,2)=gxcart(j,2)+ & gloc_sc(3,0,icg)*dtauangle(j,3,3,3) - if ((itype(3,1).ne.10).and.(nres.ge.3).and.(itype(3,1).ne.ntyp1)) & + if ((itype(3,1).ne.10).and.(nres.ge.3).and.(itype(3,molnum(3)).ne.ntyp1_molec(3))) & then gxcart(j,2)=gxcart(j,2)-gloc_sc(3,1,icg)*dtauangle(j,3,1,4) !c the - above is due to different vector direction gcart(j,2)=gcart(j,2)+gloc_sc(3,1,icg)*dtauangle(j,3,2,4) endif if (nres.gt.3) then +! if ((itype(1,1).ne.10).and.& +! ((itype(1,molnum(1)).ne.ntyp1_molec(molnum(1))))) & gxcart(j,2)=gxcart(j,2)-gloc_sc(1,1,icg)*dtauangle(j,1,1,4) !c the - above is due to different vector direction gcart(j,2)=gcart(j,2)+gloc_sc(1,1,icg)*dtauangle(j,1,2,4) @@ -3321,7 +3412,8 @@ ! write(iout,*) gloc_sc(1,1,icg),dtauangle(j,1,1,4),"gx" endif endif - if ((itype(1,1).ne.10).and.(itype(1,1).ne.ntyp1)) then + if ((itype(1,1).ne.10).and.& + (itype(1,molnum(1)).ne.ntyp1_molec(molnum(1)))) then gcart(j,2)=gcart(j,2)+gloc_sc(1,0,icg)*dtauangle(j,1,3,3) ! write(iout,*) gloc_sc(1,0,icg),dtauangle(j,1,3,3) endif @@ -3341,7 +3433,8 @@ do i=3,nres-2 do j=1,3 ! write(iout,*) "before", gcart(j,i) - if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then + if ((itype(i,1).ne.10).and.& + (itype(i,molnum(i)).ne.ntyp1_molec(molnum(i)))) then gxcart(j,i)=gxcart(j,i)+gloc_sc(2,i-2,icg) & *dtauangle(j,2,3,i+1) & -gloc_sc(1,i-1,icg)*dtauangle(j,1,1,i+2) @@ -3349,28 +3442,39 @@ *dtauangle(j,1,2,i+2) ! write(iout,*) "new",j,i, ! & gcart(j,i),gloc_sc(1,i-1,icg),dtauangle(j,1,2,i+2) - if (itype(i-1,1).ne.10) then +! if (itype(i-1,1).ne.10) then + if ((itype(i-1,1).ne.10).and.& + (itype(i-1,molnum(i-1)).ne.ntyp1_molec(molnum(i-1)))) then + gxcart(j,i)=gxcart(j,i)+gloc_sc(3,i-2,icg) & *dtauangle(j,3,3,i+1) endif - if (itype(i+1,1).ne.10) then - gxcart(j,i)=gxcart(j,i)-gloc_sc(3,i-1,icg) & +! if (itype(i+1,1).ne.10) then + if ((itype(i+1,1).ne.10).and.& + (itype(i+1,molnum(i+1)).ne.ntyp1_molec(molnum(i+1)))) then + gxcart(j,i)=gxcart(j,i)-gloc_sc(3,i-1,icg) & *dtauangle(j,3,1,i+2) gcart(j,i)=gcart(j,i)+gloc_sc(3,i-1,icg) & *dtauangle(j,3,2,i+2) endif endif - if (itype(i-1,1).ne.10) then +! if (itype(i-1,1).ne.10) then + if ((itype(i-1,1).ne.10).and.& + (itype(i-1,molnum(i-1)).ne.ntyp1_molec(molnum(i-1)))) then gcart(j,i)=gcart(j,i)+gloc_sc(1,i-2,icg)* & dtauangle(j,1,3,i+1) endif - if (itype(i+1,1).ne.10) then +! if (itype(i+1,1).ne.10) then + if ((itype(i+1,1).ne.10).and.& + (itype(i+1,molnum(i+1)).ne.ntyp1_molec(molnum(i+1)))) then gcart(j,i)=gcart(j,i)+gloc_sc(2,i-1,icg)* & dtauangle(j,2,2,i+2) ! write(iout,*) "numer",i,gloc_sc(2,i-1,icg), ! & dtauangle(j,2,2,i+2) endif - if (itype(i+2,1).ne.10) then +! if (itype(i+2,1).ne.10) then + if ((itype(i+2,1).ne.10).and.& + (itype(i+2,molnum(i+2)).ne.ntyp1_molec(molnum(i+2)))) then gcart(j,i)=gcart(j,i)+gloc_sc(2,i,icg)* & dtauangle(j,2,1,i+3) endif @@ -3380,27 +3484,32 @@ ! Setting dE/ddnres-1 if(nres.ge.4) then do j=1,3 - if ((itype(nres-1,1).ne.10).and.(itype(nres-1,1).ne.ntyp1)) then + if ((itype(nres-1,1).ne.10).and.& + (itype(nres-1,molnum(nres-1)).ne.ntyp1_molec(molnum(nres-1)))) then gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(2,nres-3,icg) & *dtauangle(j,2,3,nres) ! write (iout,*) "gxcart(nres-1)", gloc_sc(2,nres-3,icg), ! & dtauangle(j,2,3,nres), gxcart(j,nres-1) - if (itype(nres-2,1).ne.10) then - gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(3,nres-3,icg) & +! if (itype(nres-2,1).ne.10) then + if ((itype(nres-2,1).ne.10).and.& + (itype(nres-2,molnum(nres-2)).ne.ntyp1_molec(molnum(nres-2)))) then + gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(3,nres-3,icg) & *dtauangle(j,3,3,nres) endif - if ((itype(nres,1).ne.10).and.(itype(nres,1).ne.ntyp1)) then + if ((itype(nres,1).ne.10).and.& + (itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres)))) then gxcart(j,nres-1)=gxcart(j,nres-1)-gloc_sc(3,nres-2,icg) & *dtauangle(j,3,1,nres+1) gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(3,nres-2,icg) & *dtauangle(j,3,2,nres+1) endif endif - if ((itype(nres-2,1).ne.10).and.(itype(nres-2,1).ne.ntyp1)) then + if ((itype(nres-2,1).ne.10).and.& + (itype(nres-2,molnum(nres-2)).ne.ntyp1_molec(molnum(nres-2)))) then gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(1,nres-3,icg)* & dtauangle(j,1,3,nres) endif - if ((itype(nres,1).ne.10).and.(itype(nres,1).ne.ntyp1)) then + if ((itype(nres,1).ne.10).and.(itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres)))) then gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(2,nres-2,icg)* & dtauangle(j,2,2,nres+1) ! write (iout,*) "gcart(nres-1)", gloc_sc(2,nres-2,icg), @@ -3410,7 +3519,7 @@ endif ! Settind dE/ddnres if ((nres.ge.3).and.(itype(nres,1).ne.10).and. & - (itype(nres,1).ne.ntyp1))then + (itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres))))then do j=1,3 gxcart(j,nres)=gxcart(j,nres)+gloc_sc(3,nres-2,icg) & *dtauangle(j,3,3,nres+1)+gloc_sc(2,nres-2,icg) & @@ -3418,6 +3527,7 @@ enddo endif ! The side-chain vector derivatives +! print *,"gcart",gcart(:,:) return end subroutine int_to_cart #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)