changes in wham names and trial in MD

[unres4.git] / source / unres / geometry.F90

diff --git a/source/unres/geometry.F90 b/source/unres/geometry.F90
index 7af9aca..96979d5 100644 (file)
--- a/source/unres/geometry.F90
+++ b/source/unres/geometry.F90
@@ -828,7 +828,7 @@
       integer :: it1,it2,it,j
 !d    print *,' CG Processor',me,' maxgen=',maxgen
       maxsi=1000
-!      write (iout,*) 'Gen_Rand_conf: nstart=',nstart
+      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)))
@@ -840,11 +840,14 @@
           fail=.true.
           do while (fail.and.nsi.le.maxsi)
             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
@@ -858,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*))') &
@@ -906,7 +910,7 @@
         if (nsi.gt.maxsi) return 1
         endif
        call locate_next_res(i)
-!        write(iout,*) "overlap,",overlap(i-1)
+        write(iout,*) "overlap,",overlap(i-1)
        if (overlap(i-1)) then
          if (nit.lt.maxnit) then
            back=.true.
@@ -958,12 +962,13 @@
       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
@@ -974,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.
@@ -988,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)
@@ -999,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
@@ -1011,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
@@ -1096,7 +1124,7 @@
       gen_theta=theta_temp
 !     print '(a)','Exiting GENTHETA.'
       else if (mnum.eq.2) then
-       gen_theta=aa0thet_nucl(1,it,1) -0.17 + ran_number(0.0d0,0.34d0)
+       gen_theta=2.0d0 + ran_number(0.0d0,0.34d0)
       else
               gen_theta=ran_number(theta_max/2.0,theta_max)
        endif
@@ -1142,6 +1170,11 @@
         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))