X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fgeometry.F90;h=c89442f3ab68cf4f4efab0bc9a1fa347151bd5a2;hb=a18045ea1a2b2658ffe57821e33d4231012a77cf;hp=76cee6a33dbeac37a7abaab79462701bf0ce0907;hpb=51dc008cf28ea283316811e66ae2cb0b3203dc75;p=unres4.git diff --git a/source/unres/geometry.F90 b/source/unres/geometry.F90 index 76cee6a..c89442f 100644 --- a/source/unres/geometry.F90 +++ b/source/unres/geometry.F90 @@ -453,7 +453,7 @@ ! 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 + if ((itype(i,1).ne.10).and.(molnum(i).ne.5)) then vbld_inv(nres+i)=1.0d0/vbld(nres+i) else vbld_inv(nres+i)=0.0d0 @@ -845,7 +845,9 @@ enddo if (nsi.gt.maxsi) return 1 endif ! it1.ne.10 + write(iout,*) "before origin_frame" call orig_frame + write(iout,*) "after origin_frame" i=4 nstart=4 else @@ -1169,7 +1171,7 @@ return endif if (nlobit.eq.0) then - al=ran_number(0.05d0,pi/6) + al=ran_number(0.05d0,pi/2) om=ran_number(-pi,pi) return endif @@ -1543,6 +1545,7 @@ dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) + print *,i,itypi,"sc_move" dsci_inv=dsc_inv(itypi) ! do iint=1,nint_gr(i) @@ -1550,6 +1553,8 @@ if (itype(j,molnum(j)).eq.ntyp1_molec(molnum(j))) cycle ind=ind+1 itypj=iabs(itype(j,molnum(j))) + print *,j,itypj,"sc_move" + dscj_inv=dsc_inv(itypj) sig0ij=sigma(itypi,itypj) chi1=chi(itypi,itypj) @@ -1640,6 +1645,7 @@ chiom1=chi1*om1 chiom2=chi2*om2 facsig=om1*chiom1+om2*chiom2-2.0D0*om1om2*chiom12 +! print *,"TUT?",om1*chiom1,facsig,om1,om2,om12 sigsq=1.0D0-facsig*faceps1_inv sigsq_om1=(chiom1-chiom12*om2)*faceps1_inv sigsq_om2=(chiom2-chiom12*om1)*faceps1_inv @@ -3309,8 +3315,9 @@ 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) +! write(iout,*) "pierwszy gcart", gcart(j,2) if ((itype(2,1).ne.10).and.& - (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)))) then + (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)).and.(molnum(2).ne.5))) 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 @@ -3319,7 +3326,8 @@ 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,1).ne.10) then +! write(iout,*) "drugi gcart", gcart(j,2) + if((itype(2,1).ne.10).and.(molnum(2).ne.5)) 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 @@ -3359,7 +3367,7 @@ +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,1).ne.10) then + if((itype(i,1).ne.10).and.(molnum(nres-1).ne.5)) 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 @@ -3377,12 +3385,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,1).ne.10) then + if((itype(nres-2,1).ne.10).and.(molnum(nres-1).ne.5)) 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,1).ne.10) then + if((itype(nres-1,1).ne.10).and.(molnum(nres-1).ne.5)) 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) @@ -3390,7 +3398,7 @@ enddo endif ! Settind dE/ddnres-1 -!#define DEBUG +#define DEBUG #ifdef DEBUG j=1 write(iout,*)"in int to carta",nres-1,gcart(j,nres-1),gloc(nres-3,icg),dphi(j,3,nres), & @@ -3409,7 +3417,7 @@ #endif !#undef DEBUG - if(itype(nres-1,1).ne.10) then + if((itype(nres-1,1).ne.10).and.(molnum(nres-1).ne.5)) 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) @@ -3427,7 +3435,8 @@ ! The side-chain vector derivatives do i=2,nres-1 if(itype(i,1).ne.10 .and. & - itype(i,molnum(i)).ne.ntyp1_molec(molnum(i))) then + itype(i,molnum(i)).ne.ntyp1_molec(molnum(i)).and.& + (molnum(i).ne.5)) 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) @@ -3436,7 +3445,7 @@ write(iout,*)"in int to cart",i, gxcart(j,i),gloc(ialph(i,1),icg),dalpha(j,3,i), & gloc(ialph(i,1)+nside,icg),domega(j,3,i) #endif -!#undef DEBUG +#undef DEBUG enddo endif enddo @@ -3446,6 +3455,8 @@ ! INTERTYP=3 SC...Ca...Ca...SC ! calculating dE/ddc1 18 continue +! write(iout,*) "przed sccor gcart", gcart(1,2) + ! do i=1,nres ! gloc(i,icg)=0.0D0 ! write (iout,*) "poczotkoawy",i,gloc_sc(1,i,icg) @@ -3453,7 +3464,7 @@ 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,molnum(1)).ne.ntyp1_molec(molnum(1)))) then + (itype(1,molnum(1)).ne.ntyp1_molec(molnum(1))).and.(molnum(1).ne.5)) then do j=1,3 !c Derviative was calculated for oposite vector of side chain therefore ! there is "-" sign before gloc_sc @@ -3462,7 +3473,7 @@ 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,molnum(2)).ne.ntyp1_molec(molnum(2)))) then + (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2))).and.(molnum(2).ne.5)) 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)* & @@ -3481,15 +3492,17 @@ ! ommited ! & +gloc_sc(intertyp,nres-2,icg)*dtheta(j,1,3) +! write(iout,*) "przed dE/ddc2 gcart", gcart(1,2) + ! Calculating the remainder of dE/ddc2 do j=1,3 if((itype(2,1).ne.10).and. & - (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)))) then + (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)).and.(molnum(2).ne.5))) then if ((itype(1,1).ne.10).and.& - ((itype(1,molnum(1)).ne.ntyp1_molec(molnum(1)))))& + ((itype(1,molnum(1)).ne.ntyp1_molec(molnum(1)))).and.(molnum(1).ne.5))& 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,molnum(3)).ne.ntyp1_molec(3))) & + if ((itype(3,1).ne.10).and.(nres.ge.3).and.(itype(3,molnum(3)).ne.ntyp1_molec(3)).and.molnum(3).ne.5) & 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 @@ -3501,33 +3514,35 @@ 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) -! write(iout,*) gloc_sc(1,1,icg),dtauangle(j,1,2,4),"gcart" +! write(iout,*) gloc_sc(1,1,icg),dtauangle(j,1,2,4),"gcart",gcart(j,2) ! write(iout,*) gloc_sc(1,1,icg),dtauangle(j,1,1,4),"gx" endif endif if ((itype(1,1).ne.10).and.& - (itype(1,molnum(1)).ne.ntyp1_molec(molnum(1)))) then + (itype(1,molnum(1)).ne.ntyp1_molec(molnum(1))).and.(molnum(1).ne.5)) 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,1).ne.10).and.(nres.ge.3)) then + if ((itype(3,1).ne.10).and.(nres.ge.3).and.(molnum(3).ne.5)) 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,1).ne.10).and.(nres.ge.4)) then + if ((itype(4,1).ne.10).and.(nres.ge.4).and.(molnum(4).ne.5)) 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,1),itype(1,1),itype(3,1), "gcart2" enddo +! write(iout,*) "po dE/ddc2 gcart", gcart(1,2) + ! 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,1).ne.10).and.& - (itype(i,molnum(i)).ne.ntyp1_molec(molnum(i)))) then + (itype(i,molnum(i)).ne.ntyp1_molec(molnum(i))).and.(molnum(i).ne.5)) 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) @@ -3537,14 +3552,15 @@ ! & 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).and.& - (itype(i-1,molnum(i-1)).ne.ntyp1_molec(molnum(i-1)))) then + (itype(i-1,molnum(i-1)).ne.ntyp1_molec(molnum(i-1))).and.(molnum(i-1).eq.5)) 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 if ((itype(i+1,1).ne.10).and.& - (itype(i+1,molnum(i+1)).ne.ntyp1_molec(molnum(i+1)))) then + (itype(i+1,molnum(i+1)).ne.ntyp1_molec(molnum(i+1))).and.& + (molnum(i+1).ne.5)) 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) & @@ -3553,13 +3569,15 @@ endif ! 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 + (itype(i-1,molnum(i-1)).ne.ntyp1_molec(molnum(i-1))).and.& + (molnum(i-1).ne.5)) 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).and.& - (itype(i+1,molnum(i+1)).ne.ntyp1_molec(molnum(i+1)))) then + (itype(i+1,molnum(i+1)).ne.ntyp1_molec(molnum(i+1)))& + .and. (molnum(i+1).ne.5)) 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), @@ -3567,7 +3585,7 @@ endif ! 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 + (itype(i+2,molnum(i+2)).ne.ntyp1_molec(molnum(i+2))).and.(molnum(i+2).ne.5)) then gcart(j,i)=gcart(j,i)+gloc_sc(2,i,icg)* & dtauangle(j,2,1,i+3) endif @@ -3578,19 +3596,19 @@ if(nres.ge.4) then do j=1,3 if ((itype(nres-1,1).ne.10).and.& - (itype(nres-1,molnum(nres-1)).ne.ntyp1_molec(molnum(nres-1)))) then + (itype(nres-1,molnum(nres-1)).ne.ntyp1_molec(molnum(nres-1))).and.(molnum(nres-1).ne.5)) 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 if ((itype(nres-2,1).ne.10).and.& - (itype(nres-2,molnum(nres-2)).ne.ntyp1_molec(molnum(nres-2)))) then + (itype(nres-2,molnum(nres-2)).ne.ntyp1_molec(molnum(nres-2))).and.(molnum(nres-2).ne.5)) 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,molnum(nres)).ne.ntyp1_molec(molnum(nres)))) then + (itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres))).and.(molnum(nres).ne.5)) 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) & @@ -3598,11 +3616,11 @@ endif endif if ((itype(nres-2,1).ne.10).and.& - (itype(nres-2,molnum(nres-2)).ne.ntyp1_molec(molnum(nres-2)))) then + (itype(nres-2,molnum(nres-2)).ne.ntyp1_molec(molnum(nres-2))).and.(molnum(nres-2).ne.5)) 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,molnum(nres)).ne.ntyp1_molec(molnum(nres)))) then + if ((itype(nres,1).ne.10).and.(itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres))).and.(molnum(nres).ne.5)) 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), @@ -3612,13 +3630,14 @@ endif ! Settind dE/ddnres if ((nres.ge.3).and.(itype(nres,1).ne.10).and. & - (itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres))))then + (itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres))).and.(molnum(nres).ne.5))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) & *dtauangle(j,2,3,nres+1) enddo endif +! write(iout,*) "final gcart",gcart(1,2) ! The side-chain vector derivatives ! print *,"gcart",gcart(:,:) return @@ -4041,8 +4060,11 @@ do i=1,nres if (molnum(i).eq.5) then c(1,i)=dmod(c(1,i),boxxsize) + if (c(1,i).lt.0) c(1,i)=c(1,i)+boxxsize c(2,i)=dmod(c(2,i),boxysize) + if (c(2,i).lt.0) c(2,i)=c(2,i)+boxysize c(3,i)=dmod(c(3,i),boxzsize) + if (c(3,i).lt.0) c(3,i)=c(3,i)+boxzsize c(1,i+nres)=dmod(c(1,i+nres),boxxsize) c(2,i+nres)=dmod(c(2,i+nres),boxysize) c(3,i+nres)=dmod(c(3,i+nres),boxzsize)