X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fgeometry.f90;h=52c3b833785307229bd2e74bc7425d4f49c7abb5;hb=8ca97b16fe25b7053f258263899ba030572cc58f;hp=780684ecb0677df16289a16a43cffbd563bd6727;hpb=35f220f409bd5d21be33a402d79da2c23d3e0c3a;p=unres4.git diff --git a/source/unres/geometry.f90 b/source/unres/geometry.f90 index 780684e..52c3b83 100644 --- a/source/unres/geometry.f90 +++ b/source/unres/geometry.f90 @@ -1,4 +1,4 @@ - module geometry + module geometry !----------------------------------------------------------------------------- use io_units use names @@ -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 + write (iout,1212) restyp(itype(i,1),1),i,dist(i-1,i),& + alfai,be,dist(nres+i,i),rad2deg*alpha(nres+i,i,nres2+2),be1 enddo 1212 format (a3,'(',i3,')',2(f10.5,2f10.2)) @@ -314,8 +321,8 @@ real(kind=8) :: alphi,omegi,theta2 real(kind=8) :: dsci,dsci_inv,sinalphi,cosalphi,cosomegi,sinomegi real(kind=8) :: xp,yp,zp,cost2,sint2,rj -! dsci=dsc(itype(i)) -! dsci_inv=dsc_inv(itype(i)) +! dsci=dsc(itype(i,1)) +! dsci_inv=dsc_inv(itype(i,1)) dsci=vbld(i+nres) dsci_inv=vbld_inv(i+nres) #ifdef OSF @@ -352,7 +359,7 @@ xx(1)= xp*cost2+yp*sint2 xx(2)=-xp*sint2+yp*cost2 xx(3)= zp -!d print '(a3,i3,3f10.5,5x,3f10.5)',restyp(itype(i)),i, +!d print '(a3,i3,3f10.5,5x,3f10.5)',restyp(itype(i,1)),i, !d & xp,yp,zp,(xx(k),k=1,3) do j=1,3 xloc(j,i)=xx(j) @@ -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,31 +428,32 @@ 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 if (i.gt.2) then if (i.le.nres) phi(i+1)=beta(i-2,i-1,i,i+1) - if ((itype(i).ne.10).and.(itype(i-1).ne.10)) then + if ((itype(i,1).ne.10).and.(itype(i-1,1).ne.10)) then tauangle(3,i+1)=beta(i+nres-1,i-1,i,i+nres) endif - if (itype(i-1).ne.10) then + if (itype(i-1,1).ne.10) then tauangle(1,i+1)=beta(i-1+nres,i-1,i,i+1) omicron(1,i)=alpha(i-2,i-1,i-1+nres) omicron(2,i)=alpha(i-1+nres,i-1,i) endif - if (itype(i).ne.10) then + if (itype(i,1).ne.10) then 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) +! print *,i,vbld(i),"vbld(i)" vbld_inv(i)=1.0d0/vbld(i) vbld(nres+i)=dist(nres+i,i) - if (itype(i).ne.10) then + if (itype(i,1).ne.10) then vbld_inv(nres+i)=1.0d0/vbld(nres+i) else vbld_inv(nres+i)=0.0d0 @@ -503,17 +512,25 @@ #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 if (lprn) then do i=2,nres - write (iout,1212) restyp(itype(i)),i,vbld(i),& + write (iout,1212) restyp(itype(i,1),1),i,vbld(i),& rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),& rad2deg*alph(i),rad2deg*omeg(i) enddo @@ -524,7 +541,7 @@ #endif return end subroutine int_from_cart1 -#ifndef WHAM_RUN +#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER) !----------------------------------------------------------------------------- ! check_sc_distr.f !----------------------------------------------------------------------------- @@ -550,11 +567,11 @@ print *,'dv=',dv do 10 it=1,1 if (it.eq.10) goto 10 - open (20,file=restyp(it)//'_distr.sdc',status='unknown') + open (20,file=restyp(it,1)//'_distr.sdc',status='unknown') call gen_side(it,90.0D0 * deg2rad,al,om,fail) close (20) goto 10 - open (20,file=restyp(it)//'_distr1.sdc',status='unknown') + open (20,file=restyp(it,1)//'_distr1.sdc',status='unknown') do i=0,90 do j=0,72 prob(j,i)=0.0D0 @@ -642,14 +659,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 +772,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) ! @@ -808,10 +829,10 @@ maxsi=100 !d write (iout,*) 'Gen_Rand_conf: nstart=',nstart if (nstart.lt.5) then - it1=iabs(itype(2)) - phi(4)=gen_phi(4,iabs(itype(2)),iabs(itype(3))) + it1=iabs(itype(2,1)) + phi(4)=gen_phi(4,iabs(itype(2,1)),iabs(itype(3,1))) ! write(iout,*)'phi(4)=',rad2deg*phi(4) - if (nstart.lt.3) theta(3)=gen_theta(iabs(itype(2)),pi,phi(4)) + if (nstart.lt.3) theta(3)=gen_theta(iabs(itype(2,1)),pi,phi(4)) ! write(iout,*)'theta(3)=',rad2deg*theta(3) if (it1.ne.10) then nsi=0 @@ -847,9 +868,9 @@ endif return 1 endif - it1=iabs(itype(i-1)) - it2=iabs(itype(i-2)) - it=iabs(itype(i)) + it1=iabs(itype(i-1,1)) + it2=iabs(itype(i-2,1)) + it=iabs(itype(i,1)) ! print *,'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,' it2=',it2, ! & ' nit=',nit,' niter=',niter,' maxgen=',maxgen phi(i+1)=gen_phi(i+1,it1,it) @@ -929,12 +950,12 @@ nres2=2*nres data redfac /0.5D0/ overlap=.false. - iti=iabs(itype(i)) + iti=iabs(itype(i,1)) if (iti.gt.ntyp) return ! Check for SC-SC overlaps. !d print *,'nnt=',nnt,' nct=',nct do j=nnt,i-1 - itj=iabs(itype(j)) + itj=iabs(itype(j,1)) if (j.lt.i-1 .or. ipot.ne.4) then rcomp=sigmaii(iti,itj) else @@ -953,13 +974,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)) + itj=iabs(itype(j,1)) !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 +990,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 +1343,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.' @@ -1386,7 +1408,7 @@ do ires=1,ioverlap_last i=ioverlap(ires) - iti=iabs(itype(i)) + iti=iabs(itype(i,1)) if (iti.ne.10) then nsi=0 fail=.true. @@ -1453,8 +1475,8 @@ ! print *,'>>overlap_sc nnt=',nnt,' nct=',nct ind=0 do i=iatsc_s,iatsc_e - itypi=iabs(itype(i)) - itypi1=iabs(itype(i+1)) + itypi=iabs(itype(i,1)) + itypi1=iabs(itype(i+1,1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -1466,7 +1488,7 @@ do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) ind=ind+1 - itypj=iabs(itype(j)) + itypj=iabs(itype(j,1)) dscj_inv=dsc_inv(itypj) sig0ij=sigma(itypi,itypj) chi1=chi(itypi,itypj) @@ -1596,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) @@ -1807,7 +1891,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,24 +2914,42 @@ endif endif do i=1,nres-1 - iti=itype(i) - if (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0) then +! 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)).and.molnum(i).eq.1) 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,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,1).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 +! if (itype(1,1).eq.21) then ! theta(3)=90.0d0*deg2rad ! phi(4)=180.0d0*deg2rad ! vbld(2)=3.8d0 ! vbld_inv(2)=1.0d0/vbld(2) ! endif -! if (itype(nres).eq.21) then +! if (itype(nres,1).eq.21) then ! theta(nres)=90.0d0*deg2rad ! phi(nres)=180.0d0*deg2rad ! vbld(nres)=3.8d0 @@ -2857,36 +2959,49 @@ 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) + iti=itype(i,1) 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)) + + if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1 .and. di.eq.0.0d0) & + di=dsc(itype(i,molnum(i))) vbld(i+nres)=di - if (itype(i).ne.10) then + if (itype(i,1).ne.10) then vbld_inv(i+nres)=1.0d0/di 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 (iti.ne.0) then if(me.eq.king.or..not.out1file)then if (lprn) & - write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),& + 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 - iti=itype(i) + iti=itype(i,1) if(me.eq.king.or..not.out1file) & - write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),& + write (iout,'(a3,i4,7f10.3)') restyp(iti,1),i,dist(i,i-1),& rad2deg*theta(i),rad2deg*phi(i) enddo endif @@ -2920,7 +3035,7 @@ enddo enddo do i=2,nres-1 - if (itype(i).ne.10) then + if (itype(i,1).ne.10) then do j=1,3 dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i)) enddo @@ -2939,8 +3054,10 @@ cosfac=dsqrt(cosfac2) sinfac2=0.5d0/(1.0d0-costtab(i+1)) sinfac=dsqrt(sinfac2) - it=itype(i) - if (it.ne.10) then + it=itype(i,1) + + 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 @@ -2979,12 +3096,13 @@ enddo if (lprn) then do i=2,nres - iti=itype(i) + iti=itype(i,1) if(me.eq.king.or..not.out1file) & - write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),& + write (iout,'(a3,i4,3f10.5)') restyp(iti,1),i,xxref(i),& yyref(i),zzref(i) enddo endif + return end subroutine sc_loc_geom !----------------------------------------------------------------------------- @@ -2995,16 +3113,18 @@ integer :: i,j,ires,nscat real(kind=8),dimension(3,20) :: sccor real(kind=8) :: sccmj +! print *,"I am in sccenter",ires,nscat 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 return end subroutine sccenter -#ifndef WHAM_RUN +#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER) !----------------------------------------------------------------------------- subroutine bond_regular use calc_data @@ -3018,8 +3138,8 @@ do i=1,nres-1 vbld(i+1)=vbl vbld_inv(i+1)=1.0d0/vbld(i+1) - vbld(i+1+nres)=dsc(itype(i+1)) - vbld_inv(i+1+nres)=dsc_inv(itype(i+1)) + vbld(i+1+nres)=dsc(itype(i+1,1)) + vbld_inv(i+1+nres)=dsc_inv(itype(i+1,1)) ! print *,vbld(i+1),vbld(i+1+nres) enddo return @@ -3121,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).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 @@ -3134,11 +3257,11 @@ do j=1,3 gcart(j,2)=gcart(j,2)+gloc(1,icg)*dphi(j,2,4)+ & gloc(nres-2,icg)*dtheta(j,2,3)+gloc(nres-1,icg)*dtheta(j,1,4) - if(itype(2).ne.10) then + if(itype(2,1).ne.10) then gcart(j,2)=gcart(j,2)+gloc(ialph(2,1),icg)*dalpha(j,2,2)+ & gloc(ialph(2,1)+nside,icg)*domega(j,2,2) endif - if(itype(3).ne.10) then + if(itype(3,1).ne.10) then gcart(j,2)=gcart(j,2)+gloc(ialph(3,1),icg)*dalpha(j,1,3)+ & gloc(ialph(3,1)+nside,icg)*domega(j,1,3) endif @@ -3152,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).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).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 @@ -3170,11 +3297,11 @@ +gloc(i-1,icg)*dphi(j,2,i+2)+ & gloc(i,icg)*dphi(j,1,i+3)+gloc(nres+i-4,icg)*dtheta(j,2,i+1)+ & gloc(nres+i-3,icg)*dtheta(j,1,i+2) - if(itype(i).ne.10) then + if(itype(i,1).ne.10) then gcart(j,i)=gcart(j,i)+gloc(ialph(i,1),icg)*dalpha(j,2,i)+ & gloc(ialph(i,1)+nside,icg)*domega(j,2,i) endif - if(itype(i+1).ne.10) then + if(itype(i+1,1).ne.10) then gcart(j,i)=gcart(j,i)+gloc(ialph(i+1,1),icg)*dalpha(j,1,i+1) & +gloc(ialph(i+1,1)+nside,icg)*domega(j,1,i+1) endif @@ -3188,12 +3315,12 @@ dphi(j,3,nres-1)+gloc(nres-3,icg)*dphi(j,2,nres) & +gloc(2*nres-6,icg)* & dtheta(j,2,nres-1)+gloc(2*nres-5,icg)*dtheta(j,1,nres) - if(itype(nres-2).ne.10) then + if(itype(nres-2,1).ne.10) then gcart(j,nres-2)=gcart(j,nres-2)+gloc(ialph(nres-2,1),icg)* & dalpha(j,2,nres-2)+gloc(ialph(nres-2,1)+nside,icg)* & domega(j,2,nres-2) endif - if(itype(nres-1).ne.10) then + if(itype(nres-1,1).ne.10) then gcart(j,nres-2)=gcart(j,nres-2)+gloc(ialph(nres-1,1),icg)* & dalpha(j,1,nres-1)+gloc(ialph(nres-1,1)+nside,icg)* & domega(j,1,nres-1) @@ -3204,7 +3331,7 @@ do j=1,3 gcart(j,nres-1)=gcart(j,nres-1)+gloc(nres-3,icg)*dphi(j,3,nres)+ & gloc(2*nres-5,icg)*dtheta(j,2,nres) - if(itype(nres-1).ne.10) then + if(itype(nres-1,1).ne.10) then gcart(j,nres-1)=gcart(j,nres-1)+gloc(ialph(nres-1,1),icg)* & dalpha(j,2,nres-1)+gloc(ialph(nres-1,1)+nside,icg)* & domega(j,2,nres-1) @@ -3212,7 +3339,8 @@ enddo ! The side-chain vector derivatives do i=2,nres-1 - if(itype(i).ne.10 .and. itype(i).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) @@ -3230,8 +3358,9 @@ ! write (iout,*) "poczotkoawy",i,gloc_sc(1,i,icg) ! enddo if (nres.lt.2) return - if ((nres.lt.3).and.(itype(1).eq.10)) return - if ((itype(1).ne.10).and.(itype(1).ne.ntyp1)) then + if ((nres.lt.3).and.(itype(1,1).eq.10)) return + 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 @@ -3239,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).ne.10).and.(itype(2).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)* & @@ -3247,7 +3377,8 @@ endif enddo endif - if ((nres.ge.3).and.(itype(3).ne.10).and.(itype(3).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) @@ -3259,16 +3390,21 @@ ! Calculating the remainder of dE/ddc2 do j=1,3 - if((itype(2).ne.10).and.(itype(2).ne.ntyp1)) then - if (itype(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).ne.10).and.(nres.ge.3).and.(itype(3).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) @@ -3276,27 +3412,29 @@ ! write(iout,*) gloc_sc(1,1,icg),dtauangle(j,1,1,4),"gx" endif endif - if ((itype(1).ne.10).and.(itype(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 - if ((itype(3).ne.10).and.(nres.ge.3)) then + if ((itype(3,1).ne.10).and.(nres.ge.3)) then gcart(j,2)=gcart(j,2)+gloc_sc(2,1,icg)*dtauangle(j,2,2,4) ! write(iout,*) gloc_sc(2,1,icg),dtauangle(j,2,2,4) endif - if ((itype(4).ne.10).and.(nres.ge.4)) then + if ((itype(4,1).ne.10).and.(nres.ge.4)) then gcart(j,2)=gcart(j,2)+gloc_sc(2,2,icg)*dtauangle(j,2,1,5) ! write(iout,*) gloc_sc(2,2,icg),dtauangle(j,2,1,5) endif -! write(iout,*) gcart(j,2),itype(2),itype(1),itype(3), "gcart2" +! write(iout,*) gcart(j,2),itype(2,1),itype(1,1),itype(3,1), "gcart2" enddo ! If there are more than five residues if(nres.ge.5) then do i=3,nres-2 do j=1,3 ! write(iout,*) "before", gcart(j,i) - if ((itype(i).ne.10).and.(itype(i).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) @@ -3304,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).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).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).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).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).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 @@ -3335,37 +3484,42 @@ ! Setting dE/ddnres-1 if(nres.ge.4) then do j=1,3 - if ((itype(nres-1).ne.10).and.(itype(nres-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).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).ne.10).and.(itype(nres).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).ne.10).and.(itype(nres-2).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).ne.10).and.(itype(nres).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), -! & dtauangle(j,2,2,nres+1), itype(nres-1),itype(nres) +! & dtauangle(j,2,2,nres+1), itype(nres-1,1),itype(nres,1) endif enddo endif ! Settind dE/ddnres - if ((nres.ge.3).and.(itype(nres).ne.10).and. & - (itype(nres).ne.ntyp1))then + if ((nres.ge.3).and.(itype(nres,1).ne.10).and. & + (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) & @@ -3373,9 +3527,10 @@ enddo endif ! The side-chain vector derivatives +! print *,"gcart",gcart(:,:) return end subroutine int_to_cart -#ifndef WHAM_RUN +#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER) !----------------------------------------------------------------------------- ! readrtns_CSA.F !----------------------------------------------------------------------------- @@ -3431,7 +3586,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 @@ -3453,7 +3608,7 @@ write (iout,100) do i=1,nres - write (iout,110) restyp(itype(i)),i,c(1,i),c(2,i),& + write (iout,110) restyp(itype(i,1),1),i,c(1,i),c(2,i),& c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i) enddo 100 format (//' alpha-carbon coordinates ',& @@ -3479,18 +3634,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 +3692,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 !-----------------------------------------------------------------------------