Adam's unres update
[unres.git] / source / unres / src-HCD-5D / gen_rand_conf.F
index ea009b6..603e6c0 100644 (file)
@@ -13,7 +13,7 @@ C Generate random conformation or chain cut and regrowth.
       logical overlap,back,fail
       integer nstart
       integer i,j,k,it,it1,it2,nit,niter,nsi,maxsi,maxnit
-      double precision gen_theta,gen_phi,dist
+      double precision gen_theta,gen_phi,dist,ran_number
 cd    print *,' CG Processor',me,' maxgen=',maxgen
       maxsi=100
 cd    write (iout,*) 'Gen_Rand_conf: nstart=',nstart
@@ -60,12 +60,16 @@ c       write(iout,*)'theta(3)=',rad2deg*theta(3)
        it1=iabs(itype(i-1))
        it2=iabs(itype(i-2))
        it=iabs(itype(i))
-c       print *,'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,' it2=',it2,
-c    &    ' nit=',nit,' niter=',niter,' maxgen=',maxgen
+        if (it.eq.ntyp1 .and. it1.eq.ntyp1) then
+          vbld(i)=ran_number(3.8d0,10.0d0)
+          vbld_inv(i)=1.0d0/vbld(i)
+        endif
+c        write (iout,*) 'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,
+c     &    ' 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)
-c         print *,'phi(',i,')=',phi(i)
+c          write(iout,*) 'phi(',i,')=',phi(i)," type",it1,it2,it
          theta(i-1)=gen_theta(it2,phi(i-1),phi(i))
          if (it2.ne.10 .and. it2.ne.ntyp1) then
             nsi=0
@@ -78,7 +82,12 @@ c         print *,'phi(',i,')=',phi(i)
           endif
          call locate_next_res(i-1)
        endif
-       theta(i)=gen_theta(it1,phi(i),phi(i+1))
+        if (it1.ne.ntyp1) then
+         theta(i)=gen_theta(it1,phi(i),phi(i+1))
+        else
+          theta(i)=ran_number(1.326d0,2.548d0)
+        endif
+c        write (iout,*) "i",i," theta",theta(i)
         if (it1.ne.10 .and. it1.ne.ntyp1) then 
         nsi=0
         fail=.true.
@@ -86,10 +95,12 @@ c         print *,'phi(',i,')=',phi(i)
           call gen_side(it1,theta(i),alph(i-1),omeg(i-1),fail)
           nsi=nsi+1
         enddo
+c        write (iout,*) "After forward SC generation:",nsi,maxsi
         if (nsi.gt.maxsi) return1
         endif
        call locate_next_res(i)
        if (overlap(i-1)) then
+c          write (iout,*) "overlap",i-1
          if (nit.lt.maxnit) then
            back=.true.
            nit=nit+1
@@ -107,6 +118,7 @@ c         print *,'phi(',i,')=',phi(i)
            endif
           endif
         else
+c          write (iout,*) "No overlap",i-1
          back=.false.
          nit=0 
          i=i+1
@@ -211,11 +223,55 @@ c--------------------------------------------------------------------------
       double precision function gen_phi(i,it1,it2)
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
+      include 'COMMON.IOUNITS'
       include "COMMON.TORCNSTR"
       include 'COMMON.GEO'
       include 'COMMON.BOUNDS'
-      if (raw_psipred .or. ndih_constr.eq.0) then
+      include 'COMMON.INTERACT'
+      double precision sumprob(3)
+      double precision pinorm
+      external pinorm
+      if (ndih_constr.eq.0) then
         gen_phi=ran_number(-pi,pi) 
+      else if (raw_psipred) then
+        if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
+     &    .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
+          ii=iconstr_dih(i)
+          sumprob(1)=vpsipred(2,ii)
+          sumprob(2)=sumprob(1)+vpsipred(3,ii)
+          sumprob(3)=sumprob(2)+vpsipred(1,ii)
+          aux=ran_number(0.0d0,sumprob(3))
+#ifdef DEBUG
+          write(iout,*)"gen_phi: residue i",i," ii",ii," vpsipred",
+     &      vpsipred(2,ii),vpsipred(3,ii),vpsipred(1,ii)," sumprob",
+     &      sumprob(1),sumprob(2),sumprob(3)
+          write (iout,*) "aux",aux
+#endif
+          if (aux.le.sumprob(1)) then
+#ifdef DEBUG
+            write (iout,*) "hel:",
+     &        (phibound(1,i)-sdihed(1,ndih_constr))*rad2deg,
+     &        (phibound(1,i)+sdihed(1,ndih_constr))*rad2deg
+#endif
+            gen_phi=ran_number(phibound(1,i)-sdihed(1,ndih_constr)
+     &        ,phibound(1,i)+sdihed(1,ndih_constr)) 
+          else if (aux.le.sumprob(2)) then
+#ifdef DEBUG
+            write (iout,*) "ext:",
+     &          (phibound(2,i)-sdihed(2,ndih_constr))*rad2deg,
+     &          (phibound(2,i)+sdihed(2,ndih_constr))*rad2deg
+#endif
+           gen_phi=pinorm(ran_number(phibound(2,i)-sdihed(2,ndih_constr)
+     &        ,phibound(2,i)+sdihed(2,ndih_constr)))
+          else
+#ifdef DEBUG
+            write (iout,*) "coil:",-180.0,180.0
+#endif
+            gen_phi=ran_number(-pi,pi) 
+          endif
+        else
+          gen_phi=ran_number(-pi,pi) 
+        endif
       else
 C 8/13/98 Generate phi using pre-defined boundaries
         gen_phi=ran_number(phibound(1,i),phibound(2,i)) 
@@ -288,7 +344,7 @@ c-------------------------------------------------------------------------
       fail=.false.
       if (the.eq.0.0D0 .or. the.eq.pi) then
 #ifdef MPI
-        write (*,'(a,i4,a,i3,a,1pe14.5)') 
+        write (iout,'(a,i4,a,i3,a,1pe14.5)') 
      & 'CG Processor:',me,' Error in GenSide: it=',it,' theta=',the
 #else
 cd        write (iout,'(a,i3,a,1pe14.5)') 
@@ -778,23 +834,24 @@ c     overlapping residues left, or false otherwise (success)
       include 'COMMON.IOUNITS'
       logical had_overlaps,fail,scfail
       integer ioverlap(maxres),ioverlap_last
+      integer maxit_corr /5000/
 
       had_overlaps=.false.
-      call overlap_sc_list(ioverlap,ioverlap_last)
+      call overlap_sc_list(ioverlap,ioverlap_last,.false.)
       if (ioverlap_last.gt.0) then
         write (iout,*) '#OVERLAPing residues ',ioverlap_last
-        write (iout,'(18i5)') (ioverlap(k),k=1,ioverlap_last)
+        write (iout,'(15i6)') (ioverlap(k),k=1,ioverlap_last)
         had_overlaps=.true.
       endif
 
       maxsi=1000
-      do k=1,1000
+      do k=1,maxit_corr
         if (ioverlap_last.eq.0) exit
 
         do ires=1,ioverlap_last 
           i=ioverlap(ires)
           iti=iabs(itype(i))
-          if (iti.ne.10) then
+          if (iti.ne.10 .and. iti.lt.ntyp1) then
             nsi=0
             fail=.true.
             do while (fail.and.nsi.le.maxsi)
@@ -805,16 +862,15 @@ c     overlapping residues left, or false otherwise (success)
             if(fail) goto 999
           endif
         enddo
-
 c        write (iout,*) "before chaincuild overlap_sc_list: dc0",dc(:,0)
 c        call chainbuild_extconf
 c        write (iout,*) "after chaincuild overlap_sc_list: dc0",dc(:,0)
-        call overlap_sc_list(ioverlap,ioverlap_last)
-        write (iout,*) 'Overlaping residues ',ioverlap_last,
-     &           (ioverlap(j),j=1,ioverlap_last)
+        call overlap_sc_list(ioverlap,ioverlap_last,.false.)
+        write (iout,*)'#Overlaping residues @iter',k,":",ioverlap_last
+        write (iout,*)'Residue list:',(ioverlap(j),j=1,ioverlap_last)
       enddo
 
-      if (k.le.1000.and.ioverlap_last.eq.0) then
+      if (k.le.maxit_corr.and.ioverlap_last.eq.0) then
         scfail=.false.
         if (had_overlaps) then
           write (iout,*) '#OVERLAPing all corrected after ',k,
@@ -823,22 +879,29 @@ c        write (iout,*) "after chaincuild overlap_sc_list: dc0",dc(:,0)
       else
         scfail=.true.
         write (iout,*) '#OVERLAPing NOT all corrected ',ioverlap_last
-        write (iout,'(20i4)') (ioverlap(j),j=1,ioverlap_last)
+        write (iout,'(15i6)') (ioverlap(j),j=1,ioverlap_last)
       endif
 
       return
 
  999  continue
-      write (iout,'(a30,i5,a12,i4)') 
+      write (iout,'(a30,i5,a12,i6)') 
      &               '#OVERLAP FAIL in gen_side after',maxsi,
      &               'iter for RES',i
       scfail=.true.
       return
       end
 
-      subroutine overlap_sc_list(ioverlap,ioverlap_last)
-      implicit real*8 (a-h,o-z)
+      subroutine overlap_sc_list(ioverlap,ioverlap_last,lprn)
+      implicit none
       include 'DIMENSIONS'
+#ifdef MPI
+      include "mpif.h"
+      include "COMMON.SETUP"
+      integer ierror
+      integer ioverlap_last_tab(0:max_fg_procs-1),
+     & ioverlap_all(maxres*max_fg_procs),displs(0:max_fg_procs-1)
+#endif
       include 'COMMON.GEO'
       include 'COMMON.LOCAL'
       include 'COMMON.IOUNITS'
@@ -847,19 +910,47 @@ c        write (iout,*) "after chaincuild overlap_sc_list: dc0",dc(:,0)
       include 'COMMON.FFIELD'
       include 'COMMON.VAR'
       include 'COMMON.CALC'
-      logical fail
+      integer ii,itypi,itypj,itypi1,ind,ikont
+      logical fail,lprn
       integer ioverlap(maxres),ioverlap_last
-      data redfac /0.5D0/
+      double precision redfac /0.5D0/
+      double precision rrij,rij_shift,sig0ij,xi,yi,zi,rcomp,sig
+      double precision dist
 
-      write (iout,*) "overlap_sc_list"
+#ifdef MPI
+      if (nfgtasks.gt.1) then
+        if (fg_rank.eq.0) 
+     &   call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
+        call MPI_Bcast(c(1,1),6*nres,MPI_DOUBLE_PRECISION,king,FG_COMM,
+     &    IERROR)
+        call MPI_Bcast(dc(1,0),6*(nres+1),MPI_DOUBLE_PRECISION,king,
+     &   FG_COMM,IERROR)
+        call MPI_Bcast(dc_norm(1,0),6*(nres+1),MPI_DOUBLE_PRECISION,
+     &   king,FG_COMM,IERROR)
+      endif
+#endif
+c      write (iout,*) "overlap_sc_list"
 c      write(iout,*) "iatsc_s",iatsc_s," iatsc_e",iatsc_e
-      write(iout,*) "nnt",nnt," nct",nct
+c      write(iout,*) "nnt",nnt," nct",nct
       ioverlap_last=0
 C Check for SC-SC overlaps and mark residues
 c      print *,'>>overlap_sc nnt=',nnt,' nct=',nct
       ind=0
+#ifdef DEBUG
+      write (iout,*) "FG proecssor",fg_rank," g_listscsc_start",
+     & g_listscsc_start," g_listscsc_end",g_listscsc_end
+      write (*,*) "FG proecssor",fg_rank," g_listscsc_start",
+     & g_listscsc_start," g_listscsc_end",g_listscsc_end
+#endif
 c      do i=iatsc_s,iatsc_e
-      do i=nnt,nct
+      do ikont=g_listscsc_start,g_listscsc_end
+c        write (*,*) "FG processor",fg_rank," loop begins ioverlap_last",
+c     &   ioverlap_last
+        i=newcontlisti(ikont)
+        j=newcontlistj(ikont)
+c      do i=nnt,nct
+c        write (*,*) "FG processor",fg_rank," loop begins ioverlap_last",
+c     &   ioverlap_last,"ikont i j",ikont,i,j
         itypi=iabs(itype(i))
         itypi1=iabs(itype(i+1))
         if (itypi.eq.ntyp1) cycle
@@ -873,7 +964,7 @@ c      do i=iatsc_s,iatsc_e
 c
 c        do iint=1,nint_gr(i)
 c          do j=istart(i,iint),iend(i,iint)
-          do j=i+1,nct
+c          do j=i+1,nct
             ind=ind+1
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
@@ -919,10 +1010,16 @@ c     &       " sig",sig," sig0ij",sig0ij
 c            if ( rij_shift.le.0.0D0 ) then
             if ( rij_shift/sig0ij.le.0.1D0 ) then
 c              write (iout,*) "overlap",i,j
-              write (iout,'(a,i5,a,i5,a,f10.5,a,3f10.5)')
-     &         'overlap SC-SC: i=',i,' j=',j,
-     &         ' dist=',dist(nres+i,nres+j),' rcomp=',
-     &         rcomp,1.0/rij,rij_shift
+              if (lprn) then
+                write (iout,'(a,i5,a,i5,a,f10.5,a,3f10.5)')
+     &           'overlap SC-SC: i=',i,' j=',j,
+     &           ' dist=',dist(nres+i,nres+j),' rcomp=',
+     &           rcomp,1.0/rij,rij_shift
+                write (*,'(a,i2,a,i5,a,i5,a,f10.5,a,3f10.5)')
+     &           'FG processor',fg_rank,' overlap SC-SC: i=',i,' j=',j,
+     &           ' dist=',dist(nres+i,nres+j),' rcomp=',
+     &           rcomp,1.0/rij,rij_shift
+              endif
               ioverlap_last=ioverlap_last+1
               ioverlap(ioverlap_last)=i         
               do k=1,ioverlap_last-1
@@ -933,9 +1030,73 @@ c              write (iout,*) "overlap",i,j
               do k=1,ioverlap_last-1
                 if (ioverlap(k).eq.j) ioverlap_last=ioverlap_last-1
               enddo 
+c              write(*,*) "FG processor",fg_rank,i,j," ioverlap_last",
+c     &        ioverlap_last," ioverlap",(ioverlap(k),k=1,ioverlap_last)
             endif
-          enddo
+c          enddo
 c        enddo
       enddo
+#ifdef MPI
+#ifdef DEBUG
+      write (iout,*) "FG Processor",fg_rank," ioverlap_last",
+     & ioverlap_last," ioverlap",(ioverlap(i),i=1,ioverlap_last)
+      write (*,*) "FG Processor",fg_rank," ioverlap_last",ioverlap_last,
+     & " ioverlap",(ioverlap(i),i=1,ioverlap_last)
+      call flush(iout)
+#endif
+      if (nfgtasks.eq.1) return
+#ifdef DEBUG
+      write (iout,*) "Before MPI_Gather"
+      call flush(iout)
+#endif
+      call MPI_Gather(ioverlap_last,1,MPI_INTEGER,ioverlap_last_tab,
+     & 1,MPI_INTEGER,king,FG_COMM,IERROR)
+#ifdef DEBUG
+      write (iout,*) "After MPI_Gather"
+      call flush(iout)
+#endif
+#ifdef DEBUG
+      if (myrank.eq.king) 
+     & write (iout,*) "FG Processor",fg_rank,"ioverlap_last_tab",
+     & (ioverlap_last_tab(i),i=0,nfgtasks-1)
+      call flush(iout)
+#endif
+      displs(0)=0
+      do i=1,nfgtasks-1
+        displs(i)=displs(i-1)+ioverlap_last_tab(i-1)
+      enddo
+      call MPI_Gatherv(ioverlap,ioverlap_last,MPI_INTEGER,
+     & ioverlap_all,ioverlap_last_tab,displs,MPI_INTEGER,king,
+     & FG_COMM,IERROR)
+#ifdef DEBUG
+      write (iout,*) "After Gatherv"
+      call flush(iout)
+#endif
+      if (fg_rank.gt.0) return
+      ioverlap_last=0
+      do i=0,nfgtasks-1
+        ioverlap_last=ioverlap_last+ioverlap_last_tab(i)
+      enddo
+#ifdef DEBUG
+      write (iout,*) "ioverlap_last",ioverlap_last," ioverlap_last",
+     & (ioverlap_all(i),i=1,ioverlap_last)
+      call flush(iout)
+#endif
+      ii=0
+      do i=1,ioverlap_last
+        ioverlap(ii+1)=ioverlap_all(i)
+        do j=ii,1,-1
+          if (ioverlap(ii+1).eq.ioverlap(j)) goto 11
+        enddo
+        ii=ii+1
+   11   continue 
+      enddo
+      ioverlap_last=ii
+#ifdef DEBUG
+      write (iout,*) "After summing: ioverlap_last",ioverlap_last,
+     & " ioverlap",(ioverlap(i),i=1,ioverlap_last)
+      call flush(iout)
+#endif
+#endif
       return
       end