Merge branch 'UCGM' of mmka.chem.univ.gda.pl:unres4 into UCGM
[unres4.git] / source / unres / geometry.F90
index 4a96b0f..96979d5 100644 (file)
       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')
           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
       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'
       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
       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*))') &
        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.
            endif
           endif
         else
+!          write(iout,*) "tu dochodze"
          back=.false.
          nit=0 
          i=i+1
       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 
 !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.
        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)
 ! 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
       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
       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'
       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
       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)
       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'
       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)') &
         fail=.true.
         return
       endif
+      if (nlobit.eq.0) then
+         al=ran_number(0.05d0,pi/6)
+         om=ran_number(-pi,pi)
+         return
+      endif
       tant=dtan(the-pipol)
       nlobit=nlob(it)
       allocate(z(3,nlobit))
         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
       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
         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
       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)