X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fgeometry.F90;h=9641e8cbfd4519f88d3cdae9b14ae80086fa0d00;hb=6c44c356a2ec8ee148fc7a0d066a4d7317298116;hp=4a96b0f5dfaa1946b00002d78d6e6e81df15dcce;hpb=73ee8197d313a2464fae16e818d879a1940e09bf;p=unres4.git diff --git a/source/unres/geometry.F90 b/source/unres/geometry.F90 index 4a96b0f..9641e8c 100644 --- a/source/unres/geometry.F90 +++ b/source/unres/geometry.F90 @@ -566,9 +566,9 @@ dV=2.0D0*5.0D0*deg2rad*deg2rad print *,'dv=',dv do 10 it=1,1 - if (it.eq.10) goto 10 + 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 @@ -811,6 +811,7 @@ subroutine gen_rand_conf(nstart,*) ! Generate random conformation or chain cut and regrowth. use mcm_data + use random, only: iran_num,ran_number ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CHAIN' @@ -826,24 +827,27 @@ integer :: i,nstart,maxsi,nsi,maxnit,nit,niter integer :: it1,it2,it,j !d print *,' CG Processor',me,' maxgen=',maxgen - maxsi=100 -!d write (iout,*) 'Gen_Rand_conf: nstart=',nstart + maxsi=1000 + 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.(it.ne.ntyp1)) then + 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 @@ -857,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*))') & @@ -871,35 +876,41 @@ it1=iabs(itype(i-1,molnum(i-1))) it2=iabs(itype(i-2,molnum(i-2))) it=iabs(itype(i,molnum(i))) -! print *,'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,' it2=',it2, -! & ' nit=',nit,' niter=',niter,' maxgen=',maxgen + if ((it.eq.ntyp1).and.(it1.eq.ntyp1)) & + vbld(i)=ran_number(30.0D0,40.0D0) +! 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) 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) if (overlap(i-1)) then if (nit.lt.maxnit) then back=.true. @@ -918,6 +929,7 @@ endif endif else +! write(iout,*) "tu dochodze" back=.false. nit=0 i=i+1 @@ -950,12 +962,15 @@ 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 rcomp=sigmaii(iti,itj) else @@ -964,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. @@ -978,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) @@ -989,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 @@ -1001,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 @@ -1033,8 +1071,8 @@ return end function gen_phi !----------------------------------------------------------------------------- - real(kind=8) function gen_theta(it,gama,gama1) - use random,only:binorm + 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' ! include 'COMMON.LOCAL' @@ -1042,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 @@ -1060,8 +1098,13 @@ else z(1)=0.0D0 z(2)=0.0D0 - endif + endif + if (it.eq.ntyp1) then + gen_theta=ran_number(theta_max/2.0,theta_max) + else if (mnum.eq.1) then + thet_pred_mean=a0thet(it) +! write(iout,*),it,thet_pred_mean,"gen_thet" do k=1,2 thet_pred_mean=thet_pred_mean+athet(k,it,1,1)*y(k) & +bthet(k,it,1,1)*z(k) @@ -1080,10 +1123,15 @@ if (theta_temp.gt.theta_max) theta_temp=theta_max gen_theta=theta_temp ! print '(a)','Exiting GENTHETA.' + 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' @@ -1103,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)') & @@ -1121,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)) @@ -1341,16 +1395,16 @@ write(iout,'(4(1pe14.4))') ((box(k,j),k=1,2),j=1,2) write(iout,'(a,4(1pe14.4))')'sumW: ',(sumW(i),i=0,nlob(it) ) endif - if (box(1,2).lt.-MaxBoxLen .or. box(2,2).gt.MaxBoxLen) then -#ifdef MPI - 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 +! if (box(1,2).lt.-MaxBoxLen .or. box(2,2).gt.MaxBoxLen) then +!#ifdef MPI +! 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.' -#endif - fail=.true. - return - endif +!#endif +! fail=.true. +! return +! endif which_lobe=ran_number(0.0D0,sumW(nlobit)) ! print '(a,1pe14.4)','which_lobe=',which_lobe do i=1,nlobit @@ -1369,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 @@ -1409,11 +1468,11 @@ do ires=1,ioverlap_last i=ioverlap(ires) iti=iabs(itype(i,1)) - if (iti.ne.10) then + if ((iti.ne.10).and.(molnum(i).ne.5).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 @@ -1476,6 +1535,8 @@ ind=0 do i=iatsc_s,iatsc_e if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))) cycle + if (molnum(i).eq.5) print *,"WTF",i,iatsc_s,iatsc_e + if (molnum(i).eq.5) cycle itypi=iabs(itype(i,molnum(i))) itypi1=iabs(itype(i+1,1)) xi=c(1,nres+i)