corrected generation of random conformation for multichain
authorAdam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Mon, 12 Aug 2019 09:24:37 +0000 (11:24 +0200)
committerAdam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Mon, 12 Aug 2019 09:24:37 +0000 (11:24 +0200)
source/unres/MD.F90
source/unres/MREMD.F90
source/unres/geometry.F90
source/unres/io.F90

index a3451ec..a03ce09 100644 (file)
       character(len=16) :: form
       integer :: IERROR,ERRCODE
 #endif
+      integer :: iranmin,itrial,itmp
 !      include 'COMMON.SETUP'
 !      include 'COMMON.CONTROL'
 !      include 'COMMON.VAR'
           if(me.eq.king.or..not.out1file) &
              write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
          endif
+         if(iranconf.ne.0) then
+!c 8/22/17 AL Loop to produce a low-energy random conformation
+          DO iranmin=1,40
+          if (overlapsc) then
+           if(me.eq.king.or..not.out1file) &
+             write (iout,*) 'Calling OVERLAP_SC'
+           call overlap_sc(fail)
+          endif !endif overlap
+
+          if (searchsc) then
+           call sc_move(2,nres-1,10,1d10,nft_sc,etot)
+           print *,'SC_move',nft_sc,etot
+           if(me.eq.king.or..not.out1file) &
+           write(iout,*) 'SC_move',nft_sc,etot
+          endif
+
+          if(dccart)then
+           print *, 'Calling MINIM_DC'
+           call minim_dc(etot,iretcode,nfun)
+          else
+           call geom_to_var(nvar,varia)
+           print *,'Calling MINIMIZE.'
+           call minimize(etot,varia,iretcode,nfun)
+           call var_to_geom(nvar,varia)
+          endif
+          if(me.eq.king.or..not.out1file) &
+            write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
+
+          if (isnan(etot) .or. etot.gt.4.0d6) then
+            write (iout,*) "Energy too large",etot, &
+             " trying another random conformation"
+            do itrial=1,100
+              itmp=1
+              call gen_rand_conf(itmp,*30)
+              goto 40
+   30         write (iout,*) 'Failed to generate random conformation', &
+               ', itrial=',itrial
+              write (*,*) 'Processor:',me, &
+               ' Failed to generate random conformation',&
+               ' itrial=',itrial
+              call intout
+#ifdef AIX
+              call flush_(iout)
+#else
+              call flush(iout)
+#endif
+            enddo
+            write (iout,'(a,i3,a)') 'Processor:',me, &
+             ' error in generating random conformation.'
+            write (*,'(a,i3,a)') 'Processor:',me, &
+             ' error in generating random conformation.'
+            call flush(iout)
+#ifdef MPI
+!            call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
+            call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
+#else
+            stop
+#endif
+   40       continue
+          else
+            goto 44
+          endif
+          ENDDO
+
+          write (iout,'(a,i3,a)') 'Processor:',me, &
+             ' failed to generate a low-energy random conformation.'
+            write (*,'(a,i3,a,f10.3)') 'Processor:',me, &
+             ' failed to generate a low-energy random conformation.',etot
+            call flush(iout)
+            call intout
+#ifdef MPI
+!            call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
+        call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
+#else
+            stop
+#endif
+   44     continue
+        endif
       endif      
       call chainbuild_cart
       call kinetic(EK)
index 69b8ae9..dd383ad 100644 (file)
                  king,CG_COMM,ierr)
 
 ! debugging
-      print *,'traj1file',me,ii_write,ntwx_cache
+!      print *,'traj1file',me,ii_write,ntwx_cache
 ! end debugging
 
 #ifdef AIX
index 1d15584..27b2d7d 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)
         close (20)
       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'
        it1=iabs(itype(i-1,molnum(i-1)))
        it2=iabs(itype(i-2,molnum(i-2)))
        it=iabs(itype(i,molnum(i)))
+        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)
          call locate_next_res(i-1)
        endif
        theta(i)=gen_theta(it1,phi(i),phi(i+1))
-        write(iout,*) "theta(i),",theta(i)
+!        write(iout,*) "theta(i),",theta(i)
         if ((it1.ne.10).and.(it1.ne.ntyp1)) then 
         nsi=0
         fail=.true.
         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.
            endif
           endif
         else
-          write(iout,*) "tu dochodze"
+!          write(iout,*) "tu dochodze"
          back=.false.
          nit=0 
          i=i+1
 ! Check for SC-SC overlaps.
 !d    print *,'nnt=',nnt,' nct=',nct
       do j=nnt,i-1
+        
         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 
       end function gen_phi
 !-----------------------------------------------------------------------------
       real(kind=8) function gen_theta(it,gama,gama1)
-      use random,only:binorm
+      use random,only:binorm,ran_number
 !      implicit real*8 (a-h,o-z)
 !      include 'DIMENSIONS'
 !      include 'COMMON.LOCAL'
       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         
       thet_pred_mean=a0thet(it)
-      write(iout,*),it,thet_pred_mean,"gen_thet"
+!      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.'
+      endif
       return
       end function gen_theta
 !-----------------------------------------------------------------------------
         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
index 1a85795..7bcc5e2 100644 (file)
           vbld_inv(i)=vblinv
         enddo
         do i=2,nres-1
-          print *, "molnum",molnum(i),itype(i,molnum(i)),nres,i 
+!          print *, "molnum",molnum(i),itype(i,molnum(i)),nres,i 
           vbld(i+nres)=dsc(iabs(itype(i,molnum(i))))
           vbld_inv(i+nres)=dsc_inv(iabs(itype(i,molnum(i))))
 !          write (iout,*) "i",i," itype",itype(i,1),
       endif
       call read_bridge
 !--------------------------------
-       print *,"tu dochodze"
+!       print *,"tu dochodze"
 ! znamy nres oraz nss można zaalokowac potrzebne tablice
       call alloc_geo_arrays
       call alloc_ener_arrays