X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fgeometry.f90;h=52c3b833785307229bd2e74bc7425d4f49c7abb5;hb=8ca97b16fe25b7053f258263899ba030572cc58f;hp=0f3febcf12fe08af7c71f773c29cd5b1522e7b77;hpb=05fce52032feb59f7e96cd897fb82ca2aa90a888;p=unres4.git diff --git a/source/unres/geometry.f90 b/source/unres/geometry.f90 index 0f3febc..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 @@ -83,6 +83,7 @@ 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). @@ -124,7 +125,7 @@ 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),& + 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)) @@ -320,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 @@ -358,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) @@ -433,15 +434,15 @@ 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 @@ -449,9 +450,10 @@ 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 @@ -528,7 +530,7 @@ 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 @@ -565,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 @@ -827,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 @@ -866,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) @@ -948,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 @@ -976,7 +978,7 @@ 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+3).lt.4.0D0*redfac) then @@ -1406,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. @@ -1473,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) @@ -1486,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) @@ -1616,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) @@ -2850,9 +2914,11 @@ endif endif do i=1,nres-1 +! if (molnum(i).ne.1) cycle !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 + 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 @@ -2865,25 +2931,25 @@ enddo !el ----- !#ifdef WHAM_RUN -! if (itype(1).eq.ntyp1) then +! 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).eq.ntyp1) then +! 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 @@ -2897,14 +2963,15 @@ +(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 @@ -2914,18 +2981,27 @@ 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 @@ -2959,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 @@ -2978,7 +3054,7 @@ cosfac=dsqrt(cosfac2) sinfac2=0.5d0/(1.0d0-costtab(i+1)) sinfac=dsqrt(sinfac2) - it=itype(i) + it=itype(i,1) if ((it.ne.10).and.(it.ne.ntyp1)) then !el if (it.ne.10) then @@ -3020,9 +3096,9 @@ 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 @@ -3037,10 +3113,12 @@ 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 @@ -3060,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 @@ -3163,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 @@ -3176,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 @@ -3194,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 @@ -3212,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 @@ -3230,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) @@ -3246,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) @@ -3254,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) @@ -3272,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 @@ -3281,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)* & @@ -3289,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) @@ -3301,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) @@ -3318,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) @@ -3346,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 @@ -3377,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) & @@ -3415,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) @@ -3495,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 ',&