X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fgeometry.F90;h=7af9aca42933eb3d377ec2a717a5cc1875e5dcd3;hb=52d6371a0781ac6214af501666a2fbffbded3a40;hp=27b2d7deb459e76e14295fecaf3ff1040c4adf4c;hpb=a280f55746ca567cfd02f1184da61a28fbe9ebc7;p=unres4.git diff --git a/source/unres/geometry.F90 b/source/unres/geometry.F90 index 27b2d7d..7af9aca 100644 --- a/source/unres/geometry.F90 +++ b/source/unres/geometry.F90 @@ -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 @@ -833,13 +833,13 @@ 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)) nsi=nsi+1 enddo if (nsi.gt.maxsi) return 1 @@ -880,26 +880,26 @@ 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 @@ -1043,7 +1043,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 +1052,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 +1073,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 +1095,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=aa0thet_nucl(1,it,1) -0.17 + 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 +1123,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)') & @@ -1384,6 +1390,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 @@ -1428,7 +1439,7 @@ 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