random generation for nucleic acids

diff --git a/source/unres/MCM_MD.F90 b/source/unres/MCM_MD.F90
index 657f082..eb07d69 100644 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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)