X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fgeometry.f90;h=951e8eb0a5f3413bb0ed781928475c7ad17c7259;hb=9bed09ef53885a28392d806dee8108c48a7cef9d;hp=780684ecb0677df16289a16a43cffbd563bd6727;hpb=35f220f409bd5d21be33a402d79da2c23d3e0c3a;p=unres4.git diff --git a/source/unres/geometry.f90 b/source/unres/geometry.f90 index 780684e..951e8eb 100644 --- a/source/unres/geometry.f90 +++ b/source/unres/geometry.f90 @@ -9,6 +9,9 @@ use energy_data implicit none !----------------------------------------------------------------------------- +! commom.bounds +! common /bounds/ +!----------------------------------------------------------------------------- ! commom.chain ! common /chain/ ! common /rotmat/ @@ -80,17 +83,21 @@ nres2=2*nres ! Set lprn=.true. for debugging lprn = .false. + print *,"I ENTER CHAINBUILD" ! ! Define the origin and orientation of the coordinate system and locate the ! first three CA's and SC(2). ! +!elwrite(iout,*)"in chainbuild" call orig_frame +!elwrite(iout,*)"after orig_frame" ! ! Build the alpha-carbon chain. ! do i=4,nres call locate_next_res(i) enddo +!elwrite(iout,*)"after locate_next_res" ! ! First and last SC must coincide with the corresponding CA. ! @@ -111,15 +118,15 @@ write (iout,'(/a)') 'Recalculated internal coordinates' do i=2,nres-1 do j=1,3 - c(j,nres2)=0.5D0*(c(j,i-1)+c(j,i+1)) !maxres2=2*maxres + c(j,nres2+2)=0.5D0*(c(j,i-1)+c(j,i+1)) !maxres2=2*maxres enddo be=0.0D0 if (i.gt.3) be=rad2deg*beta(i-3,i-2,i-1,i) - be1=rad2deg*beta(nres+i,i,nres2,i+1) + be1=rad2deg*beta(nres+i,i,nres2+2,i+1) alfai=0.0D0 if (i.gt.2) alfai=rad2deg*alpha(i-2,i-1,i) write (iout,1212) restyp(itype(i)),i,dist(i-1,i),& - alfai,be,dist(nres+i,i),rad2deg*alpha(nres+i,i,nres2),be1 + alfai,be,dist(nres+i,i),rad2deg*alpha(nres+i,i,nres2+2),be1 enddo 1212 format (a3,'(',i3,')',2(f10.5,2f10.2)) @@ -407,6 +414,7 @@ #ifdef WHAM_RUN vbld(nres+1)=0.0d0 +!write(iout,*)"geometry warring, vbld=",(vbld(i),i=1,nres+1) vbld(2*nres)=0.0d0 vbld_inv(nres+1)=0.0d0 vbld_inv(2*nres)=0.0d0 @@ -420,7 +428,7 @@ dnorm1=dist(i-1,i) dnorm2=dist(i,i+1) do j=1,3 - c(j,nres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))/dnorm1 & + c(j,nres2+2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))/dnorm1 & +(c(j,i+1)-c(j,i))/dnorm2) enddo be=0.0D0 @@ -438,8 +446,8 @@ tauangle(2,i+1)=beta(i-2,i-1,i,i+nres) endif endif - omeg(i)=beta(nres+i,i,nres2,i+1) - alph(i)=alpha(nres+i,i,nres2) + omeg(i)=beta(nres+i,i,nres2+2,i+1) + alph(i)=alpha(nres+i,i,nres2+2) theta(i+1)=alpha(i-1,i,i+1) vbld(i)=dist(i-1,i) vbld_inv(i)=1.0d0/vbld(i) @@ -503,11 +511,19 @@ #endif do i=1,nres-1 do j=1,3 +!#ifdef WHAM_RUN +#if defined(WHAM_RUN) || defined(CLUSTER) + dc(j,i)=c(j,i+1)-c(j,i) +#endif dc_norm(j,i)=dc(j,i)*vbld_inv(i+1) enddo enddo do i=2,nres-1 do j=1,3 +!#ifdef WHAM_RUN +#if defined(WHAM_RUN) || defined(CLUSTER) + dc(j,i+nres)=c(j,i+nres)-c(j,i) +#endif dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres) enddo enddo @@ -524,7 +540,7 @@ #endif return end subroutine int_from_cart1 -#ifndef WHAM_RUN +#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER) !----------------------------------------------------------------------------- ! check_sc_distr.f !----------------------------------------------------------------------------- @@ -642,14 +658,18 @@ ii=ialph(i,2) alph(ii)=x(nphi+ntheta+i) omeg(ii)=pinorm(x(nphi+ntheta+nside+i)) +!elwrite(iout,*) "alph",ii,alph +!elwrite(iout,*) "omeg",ii,omeg enddo endif do i=4,nres phi(i)=x(i-3) +!elwrite(iout,*) "phi",i,phi enddo if (n.eq.nphi) return do i=3,nres theta(i)=x(i-2+nphi) +!elwrite(iout,*) "theta",i,theta if (theta(i).eq.pi) theta(i)=0.99d0*pi x(i-2+nphi)=theta(i) enddo @@ -751,7 +771,7 @@ thetnorm=xx return end function thetnorm -#ifndef WHAM_RUN +#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER) !----------------------------------------------------------------------------- subroutine var_to_geom_restr(n,xx) ! @@ -953,13 +973,14 @@ ! SCs. iteli=itel(i) do j=1,3 - c(j,nres2+1)=0.5D0*(c(j,i)+c(j,i+1)) +! c(j,nres2+1)=0.5D0*(c(j,i)+c(j,i+1)) + c(j,nres2+3)=0.5D0*(c(j,i)+c(j,i+1)) enddo do j=nnt,i-2 itj=iabs(itype(j)) !d print *,'overlap, p-Sc: i=',i,' j=',j, !d & ' dist=',dist(nres+j,maxres2+1) - if (dist(nres+j,nres2+1).lt.4.0D0*redfac) then + if (dist(nres+j,nres2+3).lt.4.0D0*redfac) then overlap=.true. return endif @@ -968,28 +989,28 @@ ! groups. do j=1,nnt-2 do k=1,3 - c(k,nres2+1)=0.5D0*(c(k,j)+c(k,j+1)) + c(k,nres2+3)=0.5D0*(c(k,j)+c(k,j+1)) enddo !d print *,'overlap, SC-p: i=',i,' j=',j, !d & ' dist=',dist(nres+i,maxres2+1) - if (dist(nres+i,nres2+1).lt.4.0D0*redfac) then + if (dist(nres+i,nres2+3).lt.4.0D0*redfac) then overlap=.true. return endif enddo ! Check for p-p overlaps do j=1,3 - c(j,nres2+2)=0.5D0*(c(j,i)+c(j,i+1)) + c(j,nres2+4)=0.5D0*(c(j,i)+c(j,i+1)) enddo do j=nnt,i-2 itelj=itel(j) do k=1,3 - c(k,nres2+2)=0.5D0*(c(k,j)+c(k,j+1)) + c(k,nres2+4)=0.5D0*(c(k,j)+c(k,j+1)) enddo !d print *,'overlap, p-p: i=',i,' j=',j, !d & ' dist=',dist(maxres2+1,maxres2+2) if(iteli.ne.0.and.itelj.ne.0)then - if (dist(nres2+1,nres2+2).lt.rpp(iteli,itelj)*redfac) then + if (dist(nres2+3,nres2+4).lt.rpp(iteli,itelj)*redfac) then overlap=.true. return endif @@ -1321,7 +1342,7 @@ endif if (box(1,2).lt.-MaxBoxLen .or. box(2,2).gt.MaxBoxLen) then #ifdef MPI - write (iout,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.' + write (iout,'(a,i4,a,3e15.5)') 'CG Processor:',me,': bad sampling box.',box(1,2),box(2,2),radmax write (*,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.' #else ! write (iout,'(a)') 'Bad sampling box.' @@ -1807,7 +1828,7 @@ dist=dsqrt(x12*x12+y12*y12+z12*z12) return end function dist -#ifndef WHAM_RUN +#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER) !----------------------------------------------------------------------------- ! local_move.f !----------------------------------------------------------------------------- @@ -2830,16 +2851,33 @@ endif endif do i=1,nres-1 +!in wham do i=1,nres iti=itype(i) - if (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0) 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).ne.ntyp1)) then write (iout,'(a,i4)') 'Bad Cartesians for residue',i !test stop endif +!#ifndef WHAM_RUN vbld(i+1)=dist(i,i+1) vbld_inv(i+1)=1.0d0/vbld(i+1) +!#endif if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1) if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1) enddo +!el ----- +!#ifdef WHAM_RUN +! if (itype(1).eq.ntyp1) then +! do j=1,3 +! c(j,1)=c(j,2)+(c(j,3)-c(j,4)) +! enddo +! endif +! if (itype(nres).eq.ntyp1) then +! do j=1,3 +! c(j,nres)=c(j,nres-1)+(c(j,nres-2)-c(j,nres-3)) +! enddo +! endif +!#endif ! if (unres_pdb) then ! if (itype(1).eq.21) then ! theta(3)=90.0d0*deg2rad @@ -2857,11 +2895,13 @@ if (lside) then do i=2,nres-1 do j=1,3 - c(j,nres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i) & + c(j,nres2+2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i) & +(c(j,i+1)-c(j,i))*vbld_inv(i+1)) +! in wham c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1) enddo iti=itype(i) di=dist(i,nres+i) +!#ifndef WHAM_RUN ! 10/03/12 Adam: Correction for zero SC-SC bond length if (itype(i).ne.10 .and. itype(i).ne.ntyp1 .and. di.eq.0.0d0) & di=dsc(itype(i)) @@ -2871,9 +2911,10 @@ else vbld_inv(i+nres)=0.0d0 endif +!#endif if (iti.ne.10) then - alph(i)=alpha(nres+i,i,nres2) - omeg(i)=beta(nres+i,i,nres2,i+1) + alph(i)=alpha(nres+i,i,nres2+2) + omeg(i)=beta(nres+i,i,nres2+2,i+1) endif if(me.eq.king.or..not.out1file)then if (lprn) & @@ -2940,7 +2981,9 @@ sinfac2=0.5d0/(1.0d0-costtab(i+1)) sinfac=dsqrt(sinfac2) it=itype(i) - if (it.ne.10) then + + if ((it.ne.10).and.(it.ne.ntyp1)) then +!el if (it.ne.10) then ! ! Compute the axes of tghe local cartesian coordinates system; store in ! x_prime, y_prime and z_prime @@ -2985,6 +3028,7 @@ yyref(i),zzref(i) enddo endif + return end subroutine sc_loc_geom !----------------------------------------------------------------------------- @@ -3004,7 +3048,7 @@ enddo return end subroutine sccenter -#ifndef WHAM_RUN +#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER) !----------------------------------------------------------------------------- subroutine bond_regular use calc_data @@ -3375,7 +3419,7 @@ ! The side-chain vector derivatives return end subroutine int_to_cart -#ifndef WHAM_RUN +#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER) !----------------------------------------------------------------------------- ! readrtns_CSA.F !----------------------------------------------------------------------------- @@ -3431,7 +3475,7 @@ !d & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj, !d & dhpb(i),forcon(i) !d enddo - deallocate(itype_pdb) +! deallocate(itype_pdb) return end subroutine gen_dist_constr @@ -3479,18 +3523,22 @@ ! real(kind=8),dimension(:,:),allocatable :: c !(3,maxres2+2) ! real(kind=8),dimension(:,:),allocatable :: dc allocate(dc_old(3,0:nres2)) - if(.not.allocated(dc_norm2)) allocate(dc_norm2(3,0:nres2)) !(3,0:maxres2) +! if(.not.allocated(dc_norm2)) allocate(dc_norm2(3,0:nres2+2)) !(3,0:maxres2) + if(.not.allocated(dc_norm2)) then + allocate(dc_norm2(3,0:nres2+2)) !(3,0:maxres2) + dc_norm2(:,:)=0.d0 + endif +! !el if(.not.allocated(dc_norm)) !elwrite(iout,*) "jestem w alloc geo 1" - if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:nres2)) !(3,0:maxres2) + if(.not.allocated(dc_norm)) then + allocate(dc_norm(3,0:nres2+2)) !(3,0:maxres2) + dc_norm(:,:)=0.d0 + endif !elwrite(iout,*) "jestem w alloc geo 1" allocate(xloc(3,nres),xrot(3,nres)) !elwrite(iout,*) "jestem w alloc geo 1" - do i=1,nres - do j=1,3 - xloc(j,i)=0.0D0 - enddo - enddo + xloc(:,:)=0.0D0 !elwrite(iout,*) "jestem w alloc geo 1" allocate(dc_work(6*nres)) !(MAXRES6) maxres6=6*maxres ! common /rotmat/ @@ -3533,8 +3581,17 @@ if(.not.allocated(yyref)) allocate(yyref(nres)) if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres) allocate(ialph(nres,2)) !(maxres,2) + ialph(:,1)=0 + ialph(:,2)=0 allocate(ivar(4*nres2)) !(4*maxres2) +#if defined(WHAM_RUN) || defined(CLUSTER) + allocate(vbld(2*nres)) + vbld(:)=0.d0 + allocate(vbld_inv(2*nres)) + vbld_inv(:)=0.d0 +#endif + return end subroutine alloc_geo_arrays !-----------------------------------------------------------------------------