X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fgeometry.F90;h=60c888f40fb42dd8acfd653870d0fedcc366e4e4;hb=545cf9507d923cdf917d80ed079c753702c68840;hp=27b2d7deb459e76e14295fecaf3ff1040c4adf4c;hpb=d078e021535389d9b5a48b81a90e151b8978c724;p=unres4.git diff --git a/source/unres/geometry.F90 b/source/unres/geometry.F90 index 27b2d7d..60c888f 100644 --- a/source/unres/geometry.F90 +++ b/source/unres/geometry.F90 @@ -83,7 +83,7 @@ nres2=2*nres ! Set lprn=.true. for debugging lprn = .false. - print *,"I ENTER CHAINBUILD" +! print *,"I ENTER CHAINBUILD" ! ! Define the origin and orientation of the coordinate system and locate the ! first three CA's and SC(2). @@ -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).lt.4)) then vbld_inv(nres+i)=1.0d0/vbld(nres+i) else vbld_inv(nres+i)=0.0d0 @@ -568,7 +568,7 @@ do 10 it=1,1 if ((it.eq.10).or.(it.eq.ntyp1)) goto 10 open (20,file=restyp(it,1)//'_distr.sdc',status='unknown') - call gen_side(it,90.0D0 * deg2rad,al,om,fail) + call gen_side(it,90.0D0 * deg2rad,al,om,fail,1) close (20) goto 10 open (20,file=restyp(it,1)//'_distr1.sdc',status='unknown') @@ -578,7 +578,7 @@ enddo enddo do isample=1,MaxSample - call gen_side(it,90.0D0 * deg2rad,al,om,fail) + call gen_side(it,90.0D0 * deg2rad,al,om,fail,1) indal=rad2deg*al/2 indom=(rad2deg*om+180.0D0)/5 prob(indom,indal)=prob(indom,indal)+delt @@ -828,23 +828,26 @@ integer :: it1,it2,it,j !d print *,' CG Processor',me,' maxgen=',maxgen maxsi=1000 -! write (iout,*) 'Gen_Rand_conf: nstart=',nstart + write (iout,*) 'Gen_Rand_conf: nstart=',nstart,nres if (nstart.lt.5) then 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,1)),pi,phi(4)) + if (nstart.lt.3) theta(3)=gen_theta(iabs(itype(2,1)),pi,phi(4),molnum(2)) ! write(iout,*)'theta(3)=',rad2deg*theta(3) if ((it1.ne.10).and.(it1.ne.ntyp1)) then nsi=0 fail=.true. do while (fail.and.nsi.le.maxsi) - call gen_side(it1,theta(3),alph(2),omeg(2),fail) + call gen_side(it1,theta(3),alph(2),omeg(2),fail,molnum(2)) + write (iout,*) 'nsi=',nsi,maxsi nsi=nsi+1 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 @@ -858,6 +861,7 @@ niter=0 back=.false. do while (i.le.nres .and. niter.lt.maxgen) + write(iout,*) 'i=',i,'back=',back if (i.lt.nstart) then if(iprint.gt.1) then write (iout,'(/80(1h*)/2a/80(1h*))') & @@ -880,33 +884,33 @@ if (back) then phi(i)=gen_phi(i+1,it2,it1) ! print *,'phi(',i,')=',phi(i) - theta(i-1)=gen_theta(it2,phi(i-1),phi(i)) + theta(i-1)=gen_theta(it2,phi(i-1),phi(i),molnum(i)) ! print *,"theta",theta(i-1),phi(i) if ((it2.ne.10).and.(it2.ne.ntyp1)) then nsi=0 fail=.true. do while (fail.and.nsi.le.maxsi) - call gen_side(it2,theta(i-1),alph(i-2),omeg(i-2),fail) + call gen_side(it2,theta(i-1),alph(i-2),omeg(i-2),fail,molnum(i-2)) nsi=nsi+1 enddo if (nsi.gt.maxsi) return 1 endif call locate_next_res(i-1) endif - theta(i)=gen_theta(it1,phi(i),phi(i+1)) + theta(i)=gen_theta(it1,phi(i),phi(i+1),molnum(i)) ! write(iout,*) "theta(i),",theta(i) if ((it1.ne.10).and.(it1.ne.ntyp1)) then nsi=0 fail=.true. do while (fail.and.nsi.le.maxsi) - call gen_side(it1,theta(i),alph(i-1),omeg(i-1),fail) + call gen_side(it1,theta(i),alph(i-1),omeg(i-1),fail,molnum(i)) ! write(iout,*)"alpha,omeg(i-1)",alph(i-1),omeg(i-1),i,nsi,maxsi nsi=nsi+1 enddo if (nsi.gt.maxsi) return 1 endif call locate_next_res(i) -! write(iout,*) "overlap,",overlap(i-1) + write(iout,*) "overlap,",overlap(i-1) if (overlap(i-1)) then if (nit.lt.maxnit) then back=.true. @@ -958,12 +962,13 @@ nres2=2*nres data redfac /0.5D0/ overlap=.false. - iti=iabs(itype(i,1)) + iti=iabs(itype(i,molnum(i))) if (iti.gt.ntyp) return ! Check for SC-SC overlaps. !d print *,'nnt=',nnt,' nct=',nct do j=nnt,i-1 - +! print *, "molnum(j)",j,molnum(j) + if (molnum(j).eq.1) then itj=iabs(itype(j,1)) if (itj.eq.ntyp1) cycle if (j.lt.i-1 .or. ipot.ne.4) then @@ -974,11 +979,24 @@ !d print *,'j=',j if (dist(nres+i,nres+j).lt.redfac*rcomp) then overlap=.true. + ! print *,'overlap, SC-SC: i=',i,' j=',j, ! & ' dist=',dist(nres+i,nres+j),' rcomp=', ! & rcomp return endif + else if (molnum(j).eq.2) then + itj=iabs(itype(j,2)) + if (dist(nres+i,nres+j).lt.redfac*sigma_nucl(iti,itj)) then + overlap=.true. + +! print *,'overlap, SC-SC: i=',i,' j=',j, +! & ' dist=',dist(nres+i,nres+j),' rcomp=', +! & rcomp + return + endif + + endif enddo ! Check for overlaps between the added peptide group and the preceding ! SCs. @@ -988,6 +1006,7 @@ c(j,nres2+3)=0.5D0*(c(j,i)+c(j,i+1)) enddo do j=nnt,i-2 + if (molnum(j).ne.1) cycle itj=iabs(itype(j,1)) !d print *,'overlap, p-Sc: i=',i,' j=',j, !d & ' dist=',dist(nres+j,maxres2+1) @@ -999,6 +1018,7 @@ ! Check for overlaps between the added side chain and the preceding peptide ! groups. do j=1,nnt-2 + if (molnum(j).ne.1) cycle do k=1,3 c(k,nres2+3)=0.5D0*(c(k,j)+c(k,j+1)) enddo @@ -1011,21 +1031,29 @@ enddo ! Check for p-p overlaps do j=1,3 - c(j,nres2+4)=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 +! if (molnum(j).eq.1) then itelj=itel(j) do k=1,3 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 (molnum(j).eq.1) then if(iteli.ne.0.and.itelj.ne.0)then if (dist(nres2+3,nres2+4).lt.rpp(iteli,itelj)*redfac) then overlap=.true. return endif endif + else if (molnum(j).eq.2) then + if (dist(nres2+3,nres2+4).lt.3.0) then + overlap=.true. + return + endif + endif enddo return end function overlap @@ -1043,7 +1071,7 @@ return end function gen_phi !----------------------------------------------------------------------------- - real(kind=8) function gen_theta(it,gama,gama1) + real(kind=8) function gen_theta(it,gama,gama1,mnum) use random,only:binorm,ran_number ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' @@ -1052,7 +1080,7 @@ real(kind=8),dimension(2) :: y,z real(kind=8) :: theta_max,theta_min,sig,ak !el local variables - integer :: j,it,k + integer :: j,it,k,mnum real(kind=8) :: gama,gama1,thet_pred_mean,theta_temp ! print *,'gen_theta: it=',it theta_min=0.05D0*pi @@ -1073,7 +1101,8 @@ endif if (it.eq.ntyp1) then gen_theta=ran_number(theta_max/2.0,theta_max) - else + else if (mnum.eq.1) then + thet_pred_mean=a0thet(it) ! write(iout,*),it,thet_pred_mean,"gen_thet" do k=1,2 @@ -1094,11 +1123,15 @@ if (theta_temp.gt.theta_max) theta_temp=theta_max gen_theta=theta_temp ! print '(a)','Exiting GENTHETA.' - endif + else if (mnum.eq.2) then + gen_theta=2.0d0 + ran_number(0.0d0,0.34d0) + else + gen_theta=ran_number(theta_max/2.0,theta_max) + endif return end function gen_theta !----------------------------------------------------------------------------- - subroutine gen_side(it,the,al,om,fail) + subroutine gen_side(it,the,al,om,fail,mnum) use random, only:ran_number,mult_norm1 ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' @@ -1118,13 +1151,14 @@ real(kind=8) :: Big=10.0D0 logical :: lprint,fail,lcheck !el local variables - integer :: it,i,j,k,l,nlobit,ial,iom,iii,ilob + integer :: it,i,j,k,l,nlobit,ial,iom,iii,ilob,mnum real(kind=8) :: the,al,om,detApi,wart,y2,wykl,radmax real(kind=8) :: tant,zz1,W1i,radius,zk,fac,dV,sum,sum1 real(kind=8) :: which_lobe lcheck=.false. lprint=.false. fail=.false. + if (mnum.eq.1) then if (the.eq.0.0D0 .or. the.eq.pi) then #ifdef MPI write (*,'(a,i4,a,i3,a,1pe14.5)') & @@ -1136,6 +1170,11 @@ fail=.true. return endif + if (nlobit.eq.0) then + al=ran_number(0.05d0,pi/2) + om=ran_number(-pi,pi) + return + endif tant=dtan(the-pipol) nlobit=nlob(it) allocate(z(3,nlobit)) @@ -1384,6 +1423,11 @@ if (fail) return al=y(1) om=pinorm(y(2)) + else if (mnum.eq.2) then + al=0.7+ran_number(0.0d0,0.2d0) + om=ran_number(0.0d0,3.14d0) + endif + !d print *,'al=',al,' om=',om !d stop return @@ -1424,11 +1468,11 @@ do ires=1,ioverlap_last i=ioverlap(ires) iti=iabs(itype(i,1)) - if ((iti.ne.10).and.(molnum(i).ne.5).and.(iti.ne.ntyp1)) then + if ((iti.ne.10).and.(molnum(i).lt.3).and.(iti.ne.ntyp1)) then nsi=0 fail=.true. do while (fail.and.nsi.le.maxsi) - call gen_side(iti,theta(i+1),alph(i),omeg(i),fail) + call gen_side(iti,theta(i+1),alph(i),omeg(i),fail,molnum(i)) nsi=nsi+1 enddo if(fail) goto 999 @@ -1501,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) @@ -1508,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) @@ -1598,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 @@ -1708,6 +1756,11 @@ integer :: total_ints,lower_bound,upper_bound,nint integer,dimension(0:nfgtasks) :: int4proc,sint4proc !(0:max_fg_procs) integer :: i,nexcess + if (total_ints.le.0) then + lower_bound=1 + upper_bound=0 + return + endif nint=total_ints/nfgtasks do i=1,nfgtasks int4proc(i-1)=nint @@ -2911,7 +2964,7 @@ character(len=3) :: seq,res ! character*5 atom character(len=80) :: card - real(kind=8),dimension(3,20) :: sccor + real(kind=8),dimension(3,40) :: sccor integer :: i,j,iti !el rescode, logical :: lside,lprn real(kind=8) :: di,cosfac,sinfac @@ -3130,7 +3183,7 @@ ! include 'DIMENSIONS' ! include 'COMMON.CHAIN' integer :: i,j,ires,nscat - real(kind=8),dimension(3,20) :: sccor + real(kind=8),dimension(3,40) :: sccor real(kind=8) :: sccmj ! print *,"I am in sccenter",ires,nscat do j=1,3 @@ -3267,8 +3320,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).lt.3))) 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 @@ -3277,7 +3331,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).lt.3)) 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 @@ -3317,7 +3372,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).lt.3)) 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 @@ -3335,12 +3390,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).lt.3)) 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).lt.3)) 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) @@ -3349,35 +3404,35 @@ endif ! Settind dE/ddnres-1 !#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), & - gloc(2*nres-5,icg),dtheta(j,2,nres) +!#ifdef DEBUG +! j=1 +! write(iout,*)"in int to carta",nres-1,gcart(j,nres-1),gloc(nres-3,icg),dphi(j,3,nres), & +! gloc(2*nres-5,icg),dtheta(j,2,nres) -#endif +!#endif !#undef DEBUG 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) !#define DEBUG -#ifdef DEBUG - write(iout,*)"in int to cartb",nres-1,gcart(j,nres-1),gloc(nres-3,icg),dphi(j,3,nres), & - gloc(2*nres-5,icg),dtheta(j,2,nres) - -#endif +!#ifdef DEBUG +! write(iout,*)"in int to cartb",nres-1,gcart(j,nres-1),gloc(nres-3,icg),dphi(j,3,nres), & +! gloc(2*nres-5,icg),dtheta(j,2,nres) +! +!#endif !#undef DEBUG - if(itype(nres-1,1).ne.10) then + if((itype(nres-1,1).ne.10).and.(molnum(nres-1).lt.3)) 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) !#define DEBUG -#ifdef DEBUG - write(iout,*)"in int to cart2",i,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) +!#ifdef DEBUG +! write(iout,*)"in int to cart2",i,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) -#endif +!#endif !#undef DEBUG endif @@ -3385,15 +3440,16 @@ ! 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).lt.3)) 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) !#define DEBUG -#ifdef DEBUG - 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 +!#ifdef DEBUG +! 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 enddo endif @@ -3404,6 +3460,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) @@ -3411,7 +3469,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).lt.3)) then do j=1,3 !c Derviative was calculated for oposite vector of side chain therefore ! there is "-" sign before gloc_sc @@ -3420,7 +3478,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).lt.3)) 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)* & @@ -3439,15 +3497,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).lt.3))) 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).lt.3))& 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).lt.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 @@ -3459,33 +3519,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).lt.3)) 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).lt.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,1).ne.10).and.(nres.ge.4)) then + if ((itype(4,1).ne.10).and.(nres.ge.4).and.(molnum(4).lt.3)) 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).lt.3)) 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) @@ -3495,14 +3557,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).lt.3)) 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) & @@ -3511,13 +3574,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).lt.3)) 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).lt.3)) 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), @@ -3525,7 +3590,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).lt.3)) then gcart(j,i)=gcart(j,i)+gloc_sc(2,i,icg)* & dtauangle(j,2,1,i+3) endif @@ -3536,19 +3601,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).lt.3)) 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).lt.3)) 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).lt.3)) 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) & @@ -3556,11 +3621,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).lt.3)) 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).lt.3)) 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), @@ -3570,13 +3635,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).lt.3))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 @@ -3709,7 +3775,8 @@ ! common /refstruct/ if(.not.allocated(cref)) allocate(cref(3,nres2+2,maxperm)) !(3,maxres2+2,maxperm) !elwrite(iout,*) "jestem w alloc geo 2" - allocate(crefjlee(3,nres2+2)) !(3,maxres2+2) +! allocate(crefjlee(3,nres2+2)) !(3,maxres2+2) + if (.not.allocated(crefjlee)) allocate (crefjlee(3,nres2+2)) if(.not.allocated(chain_rep)) allocate(chain_rep(3,nres2+2,maxsym)) !(3,maxres2+2,maxsym) if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym) ! common /from_zscore/ in module.compare @@ -3760,252 +3827,148 @@ !----------------------------------------------------------------------------- subroutine returnbox integer :: allareout,i,j,k,nojumpval,chain_beg,mnum - integer :: chain_end,ireturnval - real*8 :: difference + integer :: chain_end,ireturnval,idum,mnumi1 + real*8 :: difference,xi,boxsize,x,xtemp,box2shift + real(kind=8),dimension(3) :: boxx + real(kind=8),dimension(3,10000) :: xorg + integer,dimension(10000) :: posdummy + !C change suggested by Ana - end j=1 chain_beg=1 -!C do i=1,nres -!C write(*,*) 'initial', i,j,c(j,i) -!C enddo + boxx(1)=boxxsize + boxx(2)=boxysize + boxx(3)=boxzsize + idum=0 +! if(me.eq.king.or..not.out1file) then +! do i=1,nres +! write(iout,'(a6,2i3,6f9.3)') 'initial', i,j,(c(j,i),j=1,3),(c(j,i+nres),j=1,3) +! enddo +! endif !C change suggested by Ana - begin allareout=1 + chain_end=0 !C change suggested by Ana -end - do i=1,nres-1 - mnum=molnum(i) - if ((itype(i,mnum).eq.ntyp1_molec(mnum))& - .and.(itype(i+1,mnum).eq.ntyp1_molec(mnum))) then - chain_end=i - if (allareout.eq.1) then - ireturnval=int(c(j,i)/boxxsize) - if (c(j,i).le.0) ireturnval=ireturnval-1 - do k=chain_beg,chain_end - c(j,k)=c(j,k)-ireturnval*boxxsize - c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize - enddo -!C Suggested by Ana - if (chain_beg.eq.1) & - dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize -!C Suggested by Ana -end - endif - chain_beg=i+1 - allareout=1 - else - if (int(c(j,i)/boxxsize).eq.0) allareout=0 - endif - enddo - if (allareout.eq.1) then - ireturnval=int(c(j,i)/boxxsize) - if (c(j,i).le.0) ireturnval=ireturnval-1 - do k=chain_beg,nres - c(j,k)=c(j,k)-ireturnval*boxxsize - c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize - enddo - endif -!C NO JUMP -!C do i=1,nres -!C write(*,*) 'befor no jump', i,j,c(j,i) -!C enddo - nojumpval=0 - do i=2,nres - mnum=molnum(i) - if (itype(i,mnum).eq.ntyp1_molec(mnum)& - .and. itype(i-1,mnum).eq.ntyp1_molec(mnum)) then - difference=abs(c(j,i-1)-c(j,i)) -!C print *,'diff', difference - if (difference.gt.boxxsize/2.0) then - if (c(j,i-1).gt.c(j,i)) then - nojumpval=1 - else - nojumpval=-1 - endif - else - nojumpval=0 - endif - endif - c(j,i)=c(j,i)+nojumpval*boxxsize - c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize + do i=1,nres + mnum=molnum(i) + if (itype(i,mnum).eq.ntyp1_molec(mnum)) then + idum=idum+1 + posdummy(idum)=i + if (i.ne.nres) then + mnumi1=molnum(i+1) + if ((itype(i+1,mnum).eq.ntyp1_molec(mnum)).or.(mnum.ne.mnumi1)) then + do j=1,3 + xorg(j,idum)=c(j,i)-c(j,i-1) enddo - nojumpval=0 - do i=2,nres - mnum=molnum(i) - if (itype(i,mnum).eq.ntyp1_molec(mnum) .and. itype(i-1,mnum).eq.ntyp1_molec(mnum)) then - difference=abs(c(j,i-1)-c(j,i)) - if (difference.gt.boxxsize/2.0) then - if (c(j,i-1).gt.c(j,i)) then - nojumpval=1 - else - nojumpval=-1 - endif - else - nojumpval=0 - endif - endif - c(j,i)=c(j,i)+nojumpval*boxxsize - c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize + else + do j=1,3 + xorg(j,idum)=c(j,i)-c(j,i+1) enddo - -!C do i=1,nres -!C write(*,*) 'after no jump', i,j,c(j,i) -!C enddo - -!C NOW Y dimension -!C suggesed by Ana begins - allareout=1 - j=2 - chain_beg=1 - do i=1,nres-1 - mnum=molnum(i) - if ((itype(i,mnum).eq.ntyp1_molec(mnum))& - .and.(itype(i+1,mnum).eq.ntyp1_molec(mnum))) then - chain_end=i - if (allareout.eq.1) then - ireturnval=int(c(j,i)/boxysize) - if (c(j,i).le.0) ireturnval=ireturnval-1 - do k=chain_beg,chain_end - c(j,k)=c(j,k)-ireturnval*boxysize - c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize - enddo -!C Suggested by Ana - if (chain_beg.eq.1) & - dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize -!C Suggested by Ana -end - endif - chain_beg=i+1 - allareout=1 + endif else - if (int(c(j,i)/boxysize).eq.0) allareout=0 + do j=1,3 + xorg(j,idum)=c(j,i)-c(j,i-1) + enddo + endif endif enddo - if (allareout.eq.1) then - ireturnval=int(c(j,i)/boxysize) - if (c(j,i).le.0) ireturnval=ireturnval-1 - do k=chain_beg,nres - c(j,k)=c(j,k)-ireturnval*boxysize - c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize - enddo + 12 continue + do i=1,nres + mnum=molnum(i) + if (molnum(i).ge.1) then + if (i.le.3) then + k=2 + else + if (itype(i,mnum).ne.ntyp1_molec(mnum)) then + k=k+1 endif - nojumpval=0 - do i=2,nres - mnum=molnum(i) - if (itype(i,mnum).eq.ntyp1_molec(mnum)& - .and. itype(i-1,mnum).eq.ntyp1_molec(mnum)) then - difference=abs(c(j,i-1)-c(j,i)) - if (difference.gt.boxysize/2.0) then - if (c(j,i-1).gt.c(j,i)) then - nojumpval=1 - else - nojumpval=-1 - endif - else - nojumpval=0 - endif - endif - c(j,i)=c(j,i)+nojumpval*boxysize - c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize + endif +! print *,"tu2",i,k + if (itype(k,mnum).eq.ntyp1_molec(mnum)) k=k+1 + if (itype(k,mnum).eq.ntyp1_molec(mnum)) k=k+1 +! print *,"tu2",i,k + + do j=1,3 + c(j,i)=dmod(c(j,i),boxx(j)) +! if (c(1,i).lt.0) c(1,i)=c(1,i)+boxxsize + x=c(j,i)-c(j,k) +! print *,"return box,before wrap",c(1,i) + boxsize=boxx(j) + xtemp=dmod(x,boxsize) + if (dabs(xtemp-boxsize).lt.dabs(xtemp)) then + box2shift=xtemp-boxsize + else if (dabs(xtemp+boxsize).lt.dabs(xtemp)) then + box2shift=xtemp+boxsize + else + box2shift=xtemp + endif + xi=box2shift +! print *,c(1,2),xi,xtemp,box2shift + c(j,i)=c(j,k)+xi enddo - nojumpval=0 - do i=2,nres - mnum=molnum(i) - if (itype(i,mnum).eq.ntyp1_molec(mnum)& - .and. itype(i-1,mnum).eq.ntyp1) then - difference=abs(c(j,i-1)-c(j,i)) - if (difference.gt.boxysize/2.0) then - if (c(j,i-1).gt.c(j,i)) then - nojumpval=1 - else - nojumpval=-1 - endif - else - nojumpval=0 - endif - endif - c(j,i)=c(j,i)+nojumpval*boxysize - c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize + do j=1,3 + c(j,i+nres)=dmod(c(j,i+nres),boxx(j)) +! if (c(1,i).lt.0) c(1,i)=c(1,i)+boxxsize + x=c(j,i+nres)-c(j,i) +! print *,"return box,before wrap",c(1,i) + boxsize=boxx(j) + xtemp=dmod(x,boxsize) + if (dabs(xtemp-boxsize).lt.dabs(xtemp)) then + box2shift=xtemp-boxsize + else if (dabs(xtemp+boxsize).lt.dabs(xtemp)) then + box2shift=xtemp+boxsize + else + box2shift=xtemp + endif + xi=box2shift +! print *,c(1,2),xi,xtemp,box2shift + c(j,i+nres)=c(j,i)+xi enddo -!C Now Z dimension -!C Suggested by Ana -begins - allareout=1 -!C Suggested by Ana -ends - j=3 - chain_beg=1 - do i=1,nres-1 - mnum=molnum(i) - if ((itype(i,mnum).eq.ntyp1_molec(mnum))& - .and.(itype(i+1,mnum).eq.ntyp1_molec(mnum))) then - chain_end=i - if (allareout.eq.1) then - ireturnval=int(c(j,i)/boxysize) - if (c(j,i).le.0) ireturnval=ireturnval-1 - do k=chain_beg,chain_end - c(j,k)=c(j,k)-ireturnval*boxzsize - c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize - enddo -!C Suggested by Ana - if (chain_beg.eq.1) dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize -!C Suggested by Ana -end - endif - chain_beg=i+1 - allareout=1 - else - if (int(c(j,i)/boxzsize).eq.0) allareout=0 endif enddo - if (allareout.eq.1) then - ireturnval=int(c(j,i)/boxzsize) - if (c(j,i).le.0) ireturnval=ireturnval-1 - do k=chain_beg,nres - c(j,k)=c(j,k)-ireturnval*boxzsize - c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize + do i=1,idum + k=posdummy(i) + mnum=molnum(k) + if(me.eq.king.or..not.out1file) then + write(iout,*),"posdummy",i,k,(xorg(j,i),j=1,3) + endif +! do j=1,3 +! if (dabs(xorg(j,i)).gt.boxx(j))then +! x=xorg(j,i) +! boxsize=boxx(j) +! xtemp=dmod(x,boxsize) +! if (dabs(-boxsize).lt.dabs(xtemp)) then +! box2shift=xtemp-boxsize +! else if (dabs(xtemp+boxsize).lt.dabs(xtemp)) then +! box2shift=xtemp+boxsize +!! else +! box2shift=xtemp +! endif +! xorg(j,i)=box2shift +! endif +! enddo + if (k.eq.nres) then + do j=1,3 + c(j,k)=c(j,k-1)+xorg(j,i) + enddo + else + mnumi1=molnum(k+1) + if ((itype(k+1,mnum).eq.ntyp1_molec(mnum)).or.(mnum.ne.mnumi1)) then + do j=1,3 + c(j,k)=c(j,k-1)+xorg(j,i) + enddo + else + do j=1,3 + c(j,k)=c(j,k+1)+xorg(j,i) enddo - endif - nojumpval=0 - do i=2,nres - mnum=molnum(i) - if (itype(i,mnum).eq.ntyp1_molec(mnum) .and. itype(i-1,mnum).eq.ntyp1_molec(mnum)) then - difference=abs(c(j,i-1)-c(j,i)) - if (difference.gt.(boxzsize/2.0)) then - if (c(j,i-1).gt.c(j,i)) then - nojumpval=1 - else - nojumpval=-1 - endif - else - nojumpval=0 - endif - endif - c(j,i)=c(j,i)+nojumpval*boxzsize - c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize - enddo - nojumpval=0 - do i=2,nres - mnum=molnum(i) - if (itype(i,mnum).eq.ntyp1_molec(mnum) & - .and. itype(i-1,mnum).eq.ntyp1_molec(mnum)) then - difference=abs(c(j,i-1)-c(j,i)) - if (difference.gt.boxzsize/2.0) then - if (c(j,i-1).gt.c(j,i)) then - nojumpval=1 - else - nojumpval=-1 - endif - else - nojumpval=0 - endif - endif - c(j,i)=c(j,i)+nojumpval*boxzsize - c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize - enddo - do i=1,nres - if (molnum(i).eq.5) then - c(1,i)=dmod(c(1,i),boxxsize) - c(2,i)=dmod(c(2,i),boxysize) - c(3,i)=dmod(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) + endif endif enddo +! if(me.eq.king.or..not.out1file) then +! do i=1,nres +! write(iout,'(a6,2i3,6f9.3)') 'final', i,j,(c(j,i),j=1,3),(c(j,i+nres),j=1,3) +! enddo +! endif return end subroutine returnbox !-------------------------------------------------------------------------------------------------------