From 52d6371a0781ac6214af501666a2fbffbded3a40 Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Thu, 31 Oct 2019 09:57:08 +0100 Subject: [PATCH 1/1] random generation for nucleic acids --- source/unres/MCM_MD.F90 | 4 ++-- source/unres/compare.F90 | 4 ++-- source/unres/geometry.F90 | 41 ++++++++++++++++++++++++++--------------- source/unres/io.F90 | 4 ++-- source/unres/minim.F90 | 4 ++-- 5 files changed, 34 insertions(+), 23 deletions(-) diff --git a/source/unres/MCM_MD.F90 b/source/unres/MCM_MD.F90 index 657f082..eb07d69 100644 --- a/source/unres/MCM_MD.F90 +++ b/source/unres/MCM_MD.F90 @@ -3063,7 +3063,7 @@ do i=1,nside_move ii=ind_side(i) iti=itype(ii,1) - call gen_side(iti,theta(ii+1),alph(ii),omeg(ii),fail) + call gen_side(iti,theta(ii+1),alph(ii),omeg(ii),fail,1) if (fail) then isctry=isctry+1 if (isctry.gt.1000) then @@ -3084,7 +3084,7 @@ ! THETA move 400 end_select=iran_num(3,nres) theta_new=gen_theta(itype(end_select,1),phi(end_select),& - phi(end_select+1)) + phi(end_select+1),1) if (print_mc.gt.1) write (iout,'(a,i3,a,f8.3,a,f8.3)') & 'Theta ',end_select,' moved from',theta(end_select)*rad2deg,& ' to ',theta_new*rad2deg diff --git a/source/unres/compare.F90 b/source/unres/compare.F90 index eb8ea68..45ed3e5 100644 --- a/source/unres/compare.F90 +++ b/source/unres/compare.F90 @@ -4196,7 +4196,7 @@ do i=nnt,nct if (itype(i,1).ne.10) then !d print *,'i=',i,' itype=',itype(i,1),' theta=',theta(i+1) - call gen_side(itype(i,1),theta(i+1),alph(i),omeg(i),fail) + call gen_side(itype(i,1),theta(i+1),alph(i),omeg(i),fail,1) endif enddo call chainbuild @@ -4212,7 +4212,7 @@ alph0=alph(ind_sc) omeg0=omeg(ind_sc) call gen_side(itype(ind_sc,1),theta(ind_sc+1),alph(ind_sc),& - omeg(ind_sc),fail) + omeg(ind_sc),fail,1) call chainbuild call etotal(energia) !d write (iout,'(a,i5,a,i4,2(a,f8.3),2(a,1pe14.5))') 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 diff --git a/source/unres/io.F90 b/source/unres/io.F90 index ffce3d7..2a6147e 100644 --- a/source/unres/io.F90 +++ b/source/unres/io.F90 @@ -225,7 +225,7 @@ fail=.true. ii=0 do while (fail .and. ii .le. maxsi) - call gen_side(itype(i,1),theta(i+1),alph(i),omeg(i),fail) + call gen_side(itype(i,1),theta(i+1),alph(i),omeg(i),fail,molnum(i)) ii = ii+1 enddo endif @@ -1004,7 +1004,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) write(iout,*)'Adding sidechain failed for res ',& diff --git a/source/unres/minim.F90 b/source/unres/minim.F90 index 4d3115c..73b6fd2 100644 --- a/source/unres/minim.F90 +++ b/source/unres/minim.F90 @@ -4412,8 +4412,8 @@ n_try=0 do while (n_try.lt.n_maxtry .and. orig_e-cur_e.lt.e_drop) ! Move the selected residue (don't worry if it fails) - call gen_side(iabs(itype(res_pick,1)),theta(res_pick+1),& - alph(res_pick),omeg(res_pick),fail) + call gen_side(iabs(itype(res_pick,molnum(res_pick))),theta(res_pick+1),& + alph(res_pick),omeg(res_pick),fail,molnum(res_pick)) ! Minimize the side-chains starting from the new arrangement call geom_to_var(nvar,var) -- 1.7.9.5